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Interactive  Fortran  Program  to  Calculate  Thermophysical 

Properties  of  Six  Fluids 

Ben  A.  Younglove 

Thermophysical  Properties  Division 
National  Engineering  Laboratory 
National  Bureau  of  Standards 
Boulder,  Colorado  80303 

An  interactive  FORTRAN  IV  computer  program  is  given  for  computing  thermophysical 
properties  of  argon,  ethylene,  parahydrogen,  nitrogen,  nitrogen  trifluoride,  and  oxygen. 
The  program  is  designed  for  use  with  a  computer  terminal  accessing  a  large  computer  in  an 
interactive  mode.  The  program  provides  prompting  for  selection  of  several  options 
including:  I)  choice  of  fluid,  2)  choice  of  SI  or  engineering  units,  3)  choice  of  the 
single  phase  or  liquid-vapor  phase,  and  4)  a  table  of  properties  or  a  single  value. 

Properties  are  computed  for  the  single  phase  region  from  input  of  two  of  the 
variables,  temperature,  pressure,  and  density.  Values  on  the  liquid-vapor  boundary  are 
computed  from  an  entry  of  temperature  or  pressure.  The  program  returns  values  for 
pressure,  temperature,  density,  internal  energy,  enthalpy,  entropy,  specific  heats  at 
constant  volume  and  pressure,  and  sound  velocity.  Viscosity,  thermal  conductivity,  and 
dielectric  constant  are  given  for  some  of  the  fluids.  Copies  of  the  programs  may  be 
obtained  from  Office  of  Standard  Reference  Data,  Attention:  Reference  Center,  National 
Bureau  of  Standards,  Washington, D.C.  20234. 

Key  words:  argon;  computer  programs;  density;  enthalpy;  equation  of  state;  ethylene; 
hydrogen;  nitrogen;  nitrogen  trifluoride;  oxygen;  specific  heat  at  constant  pressure; 
specific  heat  at  constant  volume. 

1.  Introduction 

This  report  presents  an  interactive  FORTRAN  IV  program  for  computing  thermophysical  properties 
of  argon,  ethylene,  parahydrogen,  nitrogen,  nitrogen  trifluoride,  and  oxygen.  It  is  a  companion 
publication  to  reference  [1].  The  properties  found  in  this  reference  are  computed  from  a  standard 
set  of  equations  and  are  presented  in  tables.  References  to  the  original  works  from  which  these 
correlations  were  abstracted  are  given  in  [1], 

This  is  an  interactive  program  which  allows  selection  of  several  options  including:  1)  a  choice 
of  fluids,  2)  SI  or  engineering  units,  3)  properties  for  the  single  phase  or  for  the  liquid-vapor 
boundary,  and  4)  a  table  of  values  or  a  single  value.  Properties  of  the  homogeneous  phase  are 
computed  at  a  single  point  from  entry  of  any  two  of  the  three  variables,  temperature,  pressure,  and 
density.  Properties  on  an  isobar  may  be  specified  as  well.  Properties  on  the  liquid-vapor  boundary 
are  computed  from  an  entry  of  temperature  or  pressure. 

2.  Computer  Routines 

This  program  is  a  modification  of  program  "Fluids  Pack"  by  R.  D.  McCarty  [2].  An  early  version 
of  this  type  of  program,  which  represents  the  thermophysical  properties  of  several  fluids,  is  program 
"GASP"  (gas  properties)  by  R.  C.  Hendricks,  et  al .  [3].  The  improvements  in  the  experimental  data 
used  in  the  correlations,  that  have  occurred  since  the  publication  of  [3],  have  been  incorporated  in 
this  report.  Also,  by  retaining  the  same  form  of  the  equations  for  all  fluids,  the  computer 
programming  is  much  simplified  and  easier  to  use.  Only  the  numerical  constants,  which  are  found  in 
Appendix  E,  are  changed  from  one  fluid  to  the  next. 

A  complete  listing  of  the  computer  program  is  in  Appendices  D,  E  and  F.  The  thermophysical 
properties  produced  by  this  system  are  the  values  given  in  reference  [1].  The  program  is  designed  to 
be  run  interactively  and  provides  instructions  for  its  use  and  detects  certain  types  of  errors  such 


as  entry  of  variables  outside  of  the  range  of  validity  for  the  equation  of  state.  The  temperature 

and  pressure  range  allowed  for  each  fluid  is  automatically  displayed  v^hen  the  fluid  is  first  called. 

Since  the  basic  equation  of  state  for  the  fluids  is  semi-empirical  in  nature,  the  program  should  not 

be  used  outside  the  temperature  and  pressure  range  displayed  for  each  fluid.  Large  errors  in 

calculatea  properties  can  be  expected  outside  the  allowed  ranges.  An  example  of  the  program  in  use 

is  presented  in  appendix  B.  Properties  for  computer  program  verification  are  given  for  each  of  the 
fluids  in  appendix  C. 

3.   Thermodynamic  and  Related  Properties 

3.1  The  Equation  of  State 

The  relation  for  computing  pressure  as  a  function  of  temperature  and  density  is  a  32-term 
modified  Benedict-Webb-Rubin  (MBWR)  equation  of  state.   Its  versatility  and  adaptability  to  efficient 
computer  technique  make  it  appropriate  for  use  with  multi property  fitting  techniques,  v/here  devia- 
tions from  the  PVT  surface  and  heat  capacity  or  sound  velocity  data  are  minimized  simultaneously, 
resulting  in  a  more  accurate  model  than  would  be  the  case  if  the  PVT  data  alone  were  used.  This 
technique  has  been  described  by  R.  D.  McCarty  [2,4].  Other  advantages  of  ttie  MBWR  are:  (1)  accurate 
representation  of  the  thennodynamic  surface  over  wide  ranges  of  temperature  and  pressure, 
(2)  adaptability  to  least  squares  fitting  methods  using  many  different  kinds  of  experimental  data, 
and  (3)  convenience  in  correlating  data  from  different  sources. 

The  mathematical  form  of  the  MBWR  is, 

P  =  pRT  +  p^(G(l)T  +  G(2)T^/^  +  G(3)  +  G(4)/T  +  G(5)/T^) 
+  p^(G(6)T  +  G(7)  +  G(8)/T  +  G(9)/T^) 
+  p^(G(10)T  +  G(ll)  +  G(12)/T)  +  p^(G(13)) 
+  p^(G(14)/T  +  G(15)/t2)  +  p^(G(16)/T) 
+  p^(G(17)/T  +  G(18)t2)  +  p^(G(19)/t2) 

+  p^(G(20)/t2  +  G(21)/t2)  exp{Yp2)  (1) 

+  p^(G(22)/T^  +  G(23)/T^)  exp(YP^) 
+  p^G(24)/t2  +  G(25)/T^)  exp(YP^) 
+  p^G(26)/T^  +  G(27)/T^)  explyp^) 
+  p^^G(28)/T^  +  b(29)/T-^)  exp(Yp^) 
+  p^^(G(30)/T^  +  G(31)/T-^  +  G(32)/T^)  exp(Yp^) 

2 

The  nonlinear  coefficient,  y,  is  defined  generally  as  y  =  -1/p  ,  and  is  held  constant  during  tne 

fitting  process  to  determine  the  linear  coefficients  G(i). 

Coefficients  for  each  of  the  fluids  are  given  in  appendix  E. 

3.2  Two  Phase  Boundaries 
3.2.1  Vapor  Pressure 

The  temperature  of  the  liquid-vapor  boundary  for  each  isobar  in  the  tables  is  computed  from  a 
vapor  pressure  equation. 


?  3         4  ^(^) 

ix\?     =  £n  P^  +  Vp(l)x  +  Vp(2)x'^  +  Vp(3)x-^  +  Vp(4)x^  +  Vp(5)x(l-x)  P  (2) 

where 

X  =  (1  -  T^/T)/(i  -  yr^)  (3) 

The  coefficients  shov/n  here  as  V  ?i.rQ.   given  in  appendices  U,  E  and  F  as  VP  since  these  appendices 
?iTQ.   listings  of  Fortran  programs  and  subroutines.  The  same  applies  for  other  subscripted 
coefficients. 

3.2.2  Vapor  Densities  at  Coexistence 

Densities  of  the  vapor  in  coexistence  with  liquid,  which  are  used  to  generate  all  the  saturated 
densities  given  in  the  tables,  are  given  by 

P  =  py  (Ptv  "  Pc^  ^^P  ^^^^  ^^^ 

and 

4  -13 

f(T)  -  A(l)  Jin  X  +  E  A(i)(l  -  x^^"^)/^)  +  X)  A(i)(l  -  x)^^"^^/^         (5) 

i=2  i=5 

where 

x  =  (T  -  TJ/(T^  -  T^)   .  (6) 

3.2.3  Liquid  Densities  at  Coexistence 

The  liquid  density  at  coexistence  is  calculated  from. 


Pc  ■"  ^Pt£  "  Pc^  ^^^     ^^^^  ^^^ 


and 


f(T)  =  A(14)  £n  x  +  Y.     A(i)(l  -  x^^"^^)/^)  +  D  A(  i )  (1  -  x^^'"^^^/-^)         (8) 
i=15  i=18 

where  x  is  defined  in  eq  6. 
3.2.4  The  Melting  Line 

The  pressures  at  melting  are  given  by, 

P  (3) 
P  =  P^(l)  +  P^(2)T  ^  (9) 

3.3  Derived  Thermodynamic  Properties 

The  properties  derived  from  the  equation  of  state  and  ideal  heat  capacity  ^vo.   entropy,  enthalpy, 
internal  energy,  specific  heat  at  constant  volume  and  at  constant  pressure,  and  sound  velocity. 
3.3.1  Entropy 

The  entropy  is  computed  from 

S(T,p)  =  S^(T^)  +  J  {c°/T}dT  -  R  £n(RTp/P^)  +  J  |r/p  -  (l/p2)^|^^  I  dp   .       (10) 
0  °  P  T 


3.3.2  Ideal  Gas  Specific  Heat 

The  ideal  gas  specific  heat,  C  is  computed  from  the  following, 

C°/R  =  G.(l)/T^  +  G.(2)/T^  +  G.(3)/T  +  G.(4)  +  G.(5)T  +  G.(6)t2 


+  G,{7)t2  +  G,(8)u^e^(e'  -  1)^ 


:ii) 


where. 


3.3.3  Reference  State 


u     =     G.(9)/T     .  (12) 


The   reference  states   S     and  H     at  T  =   298.15  K  (given   below   in  table  1)   are   from  Wagman, 
et   al .    [5],   except   for  parahydrogen  whicti   is   taken   from  Woolley,    Scott,   and  Brickwedde  [5],   since 
reference  [5]   gives   only  the   norinal    hydrogen   values. 

Table  1.      Reference  Values   for  Temperature,    Entropy   and  Enthalpy. 
The  reference  value  for  temperature   is  T     =  298.15  K 
for  all    fluids.     Values  of  S°  and  H°   are  from  NBS 
TN  270-3  [5],   except   hydrogen  [6]. 

S°  J/mol    K  H°  J/mol 

argon  154.7335  6169.5 

ethylene  219.451  10564.6 

parahydrogen  130.407  8409.8 

nitrogen  191.502  8669.0 

nitrogen  trifluoride  260.621  '                 11828.0 

oxygen  2U5.029  8680.1 

3.3.4     Enthalpy 

Tlie  enthalpy   is   computed   from 


H(T,p)      -     H°  +   (P   -  pRT)/p   +/  j^  -  I^  fe)   {   dp   +  /      C°  dT     . 
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0  (  p         p     ^     /p'-j-  T^ 


3.3.5  Internal    Energy 

The   internal    energy   is, 

E(T,p)     =     H(T,p)    -  P/p      .  (14) 

3.3.6  Specific  Heat   at   Constant  Volume 

The  specific   heat   at  constant   volume   is, 

^0       „         P    U       /a^P 


C  (T,p)     =     r  -  R  -    f     \^    {^]  }     dp     .  (15) 

^  Jq     I   p         \3T 


3.3.7     Specific  Heat   at   Constant  Pressure 

The  specific   heat   as   constant   pressure   is. 


T      /3P\^,/3P 


Cp(T,p)      =     C^(T,p)    .    j^    [A)  \J^)^\     •  (16) 


3.3.8  Sound  Velocity 

The  sound  velocity  can  be  computed  as, 

W(T.P)  '\/{^)\  .  (17) 

4.  Transport  Properties:  Viscosity  and  Thermal  Conductivity 

Viscosity  and  thermal  conductivity  are   given  for  argon,  methane,  nitrogen,  and  oxygen.  The 
theory  and  discussion  of  the  relations  used  here  are  from  the  work  of  Hanley,  McCarty  and  Haynes  [7]. 
The  functional  forins  for  viscosity  and  thermal  conductivity  are, 

n  -  n^lT)  +  n^{T)p  +  n^lp.T)  (18) 

and 

^  =  Aq(T)  +  X^(T)p  +  X2(p,T)  +  A^(p,T)   .  (19) 

The  first  terms  of  eqs  (18)  and  (19)  are   the  contributions  of  the  dilute  gas, 

n<,(T)  =  i   G^(i)T<^-^)/5  , 

and  (20) 

X„(T)  =   E   G.dJTC-')/^  . 

i  =  l 

The  second  terms  in  eqs  (18)  and  (19)  represent  the  contribution  to  the  transport  coefficients  of  the 
moderately  dense  gas. 

n^d)  =  F^(l)  +  F^(2)   F^(3)  -  £n  (T/F^(4))  ^  , 
and  (21) 

A^d)  =  F^(l)  +  F^(2)   FJ3)  -  £n  (T/F^(4))  ^  . 

The  third  terms  in  these  equations  are   the  contribution  of  the  dense  gas, 

n2(p,T)  =  exp  F(p,T)  -  exp  G(T) 


:22: 


F(p,T)  =  Ey(l)  +  E^(2)H(p)  +  Ey(3)p"-^  +  £,(4)H(p)/t2 

+  Ey(5)p°-VT^-^  +  E^(6)/T  +  E^(7)H(p)/T 


:23; 


where 


G(T)  =  E^(l)  +  E^(2)/T  (24) 


H(p)  =  p°-5  (p  -  E^(8))/E^(8)   .  (25; 


The  functional  forms  of  eqs  (22),  (23),  (24)  and  (25)  are  used  for  the  corresponding  tiiermal 
conductivity  equations  with  the  coefficients  E.  substituted  for  the  E  's. 


The  last  term  of  eq  (19)  is  the  critical  enhancement  term.  Critical  effects  can  be  very  larye, 
especially  in  the  range  (p-p  )/p  <  0.25  and  (T-T  )/T  <  0.025  and  in  fact,  as  the  critical  point  is 
approached,  this  term  may  be  larger  than  the  sum  of  the  other  terms  by  several  orders  of  magnitude. 
We  include  A  to  give  the  proper  behavior  in  the  critical  region  and  to  produce  a  smooth  transition 
from  outside  the  critical  region  to  witnin  it. 


5.  Dielectric  Constant 

The  dielectric  constant  defines  the  well  known  Clausius  Mossotti  relation, 

C^  =  ^  1/p  .  (26) 

m     e+2  ' 

The  C[n  is   represented   by  the  empirical    relation, 

C^     =     B^(l)    +  B^(2)p   +  B^(3)p2  +  B^(4)p^  +  B^(5)T  +   B^(6)P     .  (27) 

Coefficients  are   given  for  parahydrogen,  nitrogen,  and  oxygen,  but  not  for  argon,  ethylene,  or 
nitrogen  trifluoride.  The  coefficients  of  Stewart  [8]  for  parahydrogen  were  modified  to  conform  to 
eq  (27).  The  data  of  Ely  and  Straty  for  nitrogen  [9]  ^^ere  used  in  a  fit  of  eq  (27).  The  oxygen 
coefficients  are  from  Younglove  [10]. 

Although  tne  range  of  the  Clausius-Mossotti  function  for  most  liquids  is  small,  the  very  high 
precision  of  the  data  require  the  terms  of  eq  (27).  The  temperature  and  pressure  terms  are   useful  in 
fitting  t\\e   compressed  liquid  states.  The  pressure  term  allows  the  flexibility  which  is  useful  in 
fluids  with  low  compressibilities,  usually  found  at  low  temperatures  and  high  pressures  [9]. 

6.  Sample  Calculations 

In  appendix  C  we  have  provided  numerical  values  computed  from  the  program  (as  given  in 
appendices  D,  E  and  F).  Values  for  the  various  properties  are   given  for  the  saturated  vapor  and 
liquid  at  a  pressure  of  0.101325  MPa  (1  atmosphere),  and  for  the  single  phase  at  35  MPa  and  300  K. 
Failure  to  reproduce  these  numbers  usually  arises  from  using  too  few  significant  figures.  This  work 
was  done  using  a  60  bit  word.  Use  of  double  precision  calculations  may  be  necessary  on  computers 
with  smaller  word  length. 
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Appendix  A.     List  of  Symbols  and  Units. 
Primary  Thermophysical   Quantities. 

P,P^,P^  pressure,   P  at  critical    point,   P  at  triple   point,   MPa 

p  ,p  density,  p   at   critical    point,  mol/dm'^ 

Pt£'Ptv  P  °^  liquid   at   triple   point,  p   of  gas   at   triple  point 

T,T_,T^  temperature,   T  at  critical    point,   T  at  triple   point,  K 

Cp,Cp  specific   heat   at   constant   pressure,   C     of  ideal    gas,  J/mol    K 

Cy,C^  specific   heat   at   constant  volume,   C     for   ideal    gas,  J/mol   K 

E  internal    energy,   J/mol 

H,H°  enthalpy,   reference  value  of  H,  J/mol 

S,S°  entropy,    reference  value  of  S,  J/mol   K 

W  sound   velocity,  m/s 

n  viscosity,  yPa  s 

A  thermal    conductivity,   W/m  K 

e  dielectric  constant 

Other  Variables  and  Constants. 

R  the  gas   constant,  8.31434  J/mol   K 

K4.  isothermal    bulk  modulus,   (MPa)"-'- 

M|„  molecular  weight,  g/mol 

Ng  Avagodro's   number,   per  mol 

C^  CI ausius-Mossotti   function 

T  reference  temperature,   298.15   K 

Critical    Point  Values. 

P(,    (MPa)  p^   (mol/dm^)  T^   (K) 

Ar  4.9058  13.41        •      150.86 

C2H4  5.0404  7.650  282.34 

H2  1.28377  15.556  32.938 

N2  3.39908  11.21  126.26 

Nf3  4.4607  7.92  234.0 

O2  5.043  13.630  154.581 


Appendix  B.  Sample  Interactive  Session 


CALL, FLUIDS 

FOR  INFORMATION  ON  HOW  TO  USE  THIS  PROGRAM,  ENTER  "0"  OTHERWISE,  ENTER  "1' 
?  0 

WHEN  THE  PROGRAM  ASKS  FOR  A  FLUID  SELECTION,  ENTER  THE  APPROPRIATE 

NUMBER.  AN  INAPPROPRIATE  NUMBER  WILL  TERMINATE  THE  PROGRAM.  WHEN 

IT  ASKS  FOR  A  PRESSURE,  DENSITY,  AND  TEMPERATURE,  ENTER  ANY  TWO  OF 

OF  THE  THREE  AND  A  ZERO  FOR  THE  THIRD. 

THE  ORDER  MUST  BE  P,D,T,  AND  ONE  OF  THE  THREE  MUST  BE  ZERO. 

IF  ALL  THREE  ARE  ZERO  THE  PROGRAM  ASKS  FOR  A  NEW  FLUID. 

SELECT  A  FLUID  BY  ENTERING  THE  CORRESPONDING  NUMBER 
1=ARG0N,      2=ETHYLENE,  3=PARA  HYDROGEN, 

4=NITR0GEN,   5=NITR0GEN  TRIFLUORIDE,   6=0XYGEN  7=ST0P 

?  1 

THE  TEMPERATURE  RANGE  FOR  ARGON  IS  83.8  TO  400  K 
WITH  PRESSURES  TO  100  MPA 

FOR  ENGINEERING  UNITS  ENTER  "0",  FOR  METRIC  ENTER  "1" 
?  1 

FOR  SATURATION  PROPERTIES  ENTER  "0",  FOR  FLUID  ENTER  "1" 
?  1 

FOR  A  SINGLE  POINT  ENTER  "0",  FOR  A  TABLE  ENTER  "1" 
?  1 

ENTER  A  PRESSURE(MPA),  A  STARTING  TEMPERATURE(K) ,  A  FINAL  TEMPERATURE 
AND  A  TEMPERATURE  INCREMENT  AND  IN  THAT  ORDER. 
?  10  100  150  10 

P     T      DEN  E      H      S    CV    CP  SOUND  VISC  COND 

MPA   K      MOL/L  ===  J/MOL/  ==  ====  J/MOL-K  -==  M/S  PA-S  MW/M-K 

E+6 

10.000  100.000  33.744  -4271.  -3974.  59.6  20.4  44.0  800.  205.0  115.1 

10.000  110.000  32.191  -3839.  -3529.  63.8  19.5  45.2  739.  162.7  102.9 

10.000  120.000  30.525  -3396.  -3069.  67.8  18.8  46.9  672.  131.2  92.1 

10.000  130.000  28.695  -2936.  -2587.  71.7  18.2  49.5  602.  106.4  81.8 

10.000  140.000  26.609  -2447.  -2072.  75.5  17.6  54.1  526.  85.9  71.7 

10.000  150.000  24.096  -1909.  -1494.  79.5  17.3  62.3  444.  67.9  61.7 

ENTER  A  PRESSURE(MPA),  A  STARTING  TEMPERATURE(K) ,  A  FINAL  TEMPERATURE 
AND  A  TEMPERATURE  INCREMENT  AND  IN  THAT  ORDER 
?  0  0  0  0 

SELECT  A  FLUID  BY  ENTERING  THE  CORRESPONDING  NUMBER 
1=ARG0N,      2=ETHYLENE,  3=PARA  HYDROGEN, 

4=NITR0GEN,   5=NITR0GEN  TRIFLUORIDE,   6=0XYGEN  7=ST0P 

?  8 

EXIT. 


Appendix  C.     Properties   for  Computer  Program  Verification 


ARGON   PROPERTIES. 


T              DEN                  E                H              S          CV            CP  SOUND  VI SC      COND 

K              MOL/L          ===  J/MOL  ==-      ====  J/MOL-K  =--=  M/S  PA-S  MW/M-K 

E+6 

SATURATED  VAPOR  AT  0.101325  MPA. 

87.2815   .1446  1038.7  1739.4  128.63  13.22  22.76  170.7  7.3   5.9 

SATURATED  LIQUID  AT  0.101325  MPA. 

87.2815  34.8946  -4708.7  -4705.8  54.82  22.27  44.73  814.4  260.3  128.3 

SINGLE  PHASE  AT  35  MPA  AND  300  K. 

300.0000  13.6426  1830.3  4395.8  100.94  14.07  32.82  445.8  39.5  36.9 


ETHYLENE  PROPERTIES 


T      DEN       E      H      S    CV  CP    SOUND  VISC  COND 

K      MOL/L    ===  J/MOL  =  =  =  ====  J/MOL-K  ====   M./S    PA-S  MW/M-K 

E+6 

SATURATED  VAPOR  AT  0.101325  MPA. 

169.4093   .0744  4180.6  5541.6  197.66  27.13  36.89  252.3 

SATURATED  LIQUID  AT  0.101325  MPA. 

169.4093  20.2449  -8001.8  -7996.8  117.73  38.22  67.92  1310.1 

SINGLE  PHASE  AT  35  MPA  AND  300  K. 

300.0000  15.8136  -504.6  1708.7  151.74  41.68  70.79  868.6 


HYDROGEN  PROPERTIES 


T    DEN       E       H      S    CV     CP    SOUND  VISC  COND  DIEL 
K    MOL/L   ===  J/MOL  ===  --==  J/MOL-K  -===   M/S   PA-S  MW/M-K  == 

E+6 

SATURATED  VAPOR  AT  0.101325  MPA. 

20.27   .564   22^.5   382.0  60.40  13.31  24.68  354.  1.00404 

SATURATED  LIQUID  AT  0.101325  MPA. 

20.27  35.118  -518.9  -516.0  16.11  11.43  19.49  1100.  1.22981 

SINGLE  PHASE  AT  35  MPA  AND  300  K. 
300.00  11.517  5837.0  8875.8  81.58  22.23  31.12  1609.  1.07174 


NITROGEN  PROPERTIES, 


T    DEN       E       H      S    CV     CP  SOUND  VISC  COND  DIEL 

K    MOL/L   ===  J/MOL  ===  ====  J/MOL-K  ====  M/S   PA-S  MW/M-K  == 

E+6 

SATURATED  VAPOR  AT  0.101325  MPA. 

77.36   .154  1548.0  2163.5  151.54  21.59  31.45  175.   5.3   7.5  1.00217 

SATURATED  LIQUID  AT  0.101325  MPA. 

77.36  28.863  -3404.3  -3400.8  79.54  27.82  57.80  939.  151.6  133.7  1.43356 

SINGLE  PHASE  AT  35  MPA  AND  300  K. 

300.00  11.756  4623. q  7601.1  138.73  22.29  38.11  514.  28.7  48.0  1.16337 
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NITROGEN  TRIFLUORIDE  PROPERTIES. 


T      DEN       E       H      S    CV     CP  SOUND  VISC  COND 

K      MOL/L    ===  J/MOL  ===  ====  J/MOL-K  ====  M/S    PA-S  MW/M-K 

E+6 
SATURATED  VAPOR  AT  0.10132S  MPA. 

144.1082   .0879  3575.7  4728.5  227.83  29.65  39.49  144.1 

SATURATED  LIQUID  AT  0.101325  MPA. 

144.1082  21.6563  -6832.5  -6827.8  147.66  40.30  72.20  757.5 

SINGLE  PHASE  AT  35  MPA  AND  300  K. 

300.0000  14.3859  3299.9  5732.9  197.40  50.03  80.58  412.4 


OXYGEN  PROPERTIES. 


T    DEN       E      H      S    CV  CP  SOUND  VISC  COND  DIEL 

K    MOL/L   ===  J/MOL  ===  ====  J/MOL-K  ====  M/S   PA-S  MW/M-K  == 

E+6 

SATURATED  VAPOR  AT  0.101325  MPA. 

90.19   .140  1811.5  2535.6  169.58  21.72  31.34  178.   6.7   8.5  1.00166 

SATURATED  LIQUID  AT  0.101325  MPA. 

90.19  35.658  -4265.5  -4262.7  94.19  29.54  54.22  905.  194.7  151.8  1.48704 

SINGLE  PHASE  AT  35  MPA  AND  300  K. 

300.00  13.742  4233.0  6780.0  150.97  22.62  41.95  442.  34.7  46.4  1.17222 
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Appendix  D.  Program  Listing  for  Program  Fluids 


PROGRAM  FLUIDS! INPUT, OUTPUT) 

DIMENSION  G(32),VP(9) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

COMMON/CONT/IF 

COMMON/CRIT/EM,EOK,RM,TC,DC,X,PC,SIG 

C0MM0N/DIEL/BX(6),PX(6) 
1000  FORMAT(*  FOR  INFORMATION  ON  HOW  TO  USE  THIS  PROGRAM,  ENTER  "0"  * 
A  *  OTHERWISE,  ENTER  "1"*) 

PRINT  1000 

IP  =  3 

PRINT  1210 

READ  1010,1 

IF{I.EQ.O)CALL  INFO 

PRINT  1210 

1010  FORMAT(Il) 

110  PRINT  1020 

PRINT  1030 
1020  FORMAT(*  SELECT  A  FLUID  BY  ENTERING  THE  CORRESPONDING  NUMBER*) 
1030  FORMAT(*  1=ARG0N,      2=ETHYLENE,  3=PARA  HYDROGEN, 

A    *,/,*  4=NITR0GEN,   5=NITR0GEN  TRIFLUORIDE,   6-OXYGEN* 
B  ,*  7=ST0P*) 

READ  *,IF 

PRINT  1210 

GO  T0(1,2,3,4,5,6,999),IF 

1  CALL  DATA  AR 
GO  TO  120 

2  CALL  DATA  C3 
GO  TO  120 

3  CALL  DATA  PH2 
GO  TO  120 

4  CALL  DATA  N2 
GO  TO  120 

5  CALL  DATA  NF3 
GO  TO  120 

6  CALL  DATA  02 
GO  TO  120 

120  PRINT  1040 
1040  FORMAT(*  FOR  ENGINEERING  UNITS  ENTER  "0",  FOR  METRIC  ENTER  "1"*) 
READ  1010,  lU 
PRINT  1210 
PRINT  1050 
1050  FORMAT(*  FOR  SATURATION  PROPERTIES  ENTER  "0",  FOR  FLUID  ENTER  * 
A  *"1"*) 
READ  1010, IC 
PRINT  1210 
PRINT  1060 
1060  FORMAT(*  FOR  A  SINGLE  POINT  ENTER  "0",  FOR  A  TABLE  ENTER  "1"*) 
READ  1010, IV 
PRINT  1210 
160  IF(IC.EQ.0)G0  TO  240 
IFdV.EQ.DGO  TO  330 
170  IF(IU.EQ.O)GO  TO  180 
PRINT  1080 
READ  *,P,D,T 
PRINT  1210 
GO  TO  190 
180  PRINT  1070 
1070  FORMAT!*  ENTER  PRESSURE! PSI A) ,  DENSITY! LB/CU  FT),  AND  TEMPERATURE* 
A  *!F)*) 
READ  *,P,D,T 
PRINT  1210 
P=!P/14. 695949)*. 101325 
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D=D*16.01846371/EM 

IF(T.E0.0.0)G0  TO  190 

T=(T-32.)/l. 8+273. 15 
190  IF{P.LE.O.O)GO  TO  220 
1080  FORMAT(*  ENTER  PRESSURE(MPA) ,  DENSITY(MOL/L) ,  AND  TEMPERATURE(K)^ 

IF(D.LE.O.O)GO  TO  210 

IF(T.LE.0.0)G0  TO  200 

GO  TO  170 
200  IF(P.LE.0.0R.D.LE.0.0)G0  TO  110 

T=FIND  T(P,D) 

CALL  LIMITS(P,T,IL) 

IF(IL.LE.O)GO  TO  170 

GO  TO  230 
210  IF(T.LE.n.0R.P.LE.0)G0  TO  110 

CALL  LIMITS(P,T,IL) 

IF(IL.LE.0)G0  TO  170 

D=FIND  D(P,T) 

GO  TO  230 
220  IF(D.LE.O.0.OR.T.LE.O)G0  TO  110 

P=FIND  P(D,T) 

CALL  LIMITS(P,T,IL) 

IF(IL.LE.O)GO  TO  170 
230  CALL  REPRO(P,D,T,IU,IV,IC,IP,TF,DELT) 

GO  TO  170 
240  PRINT  1090 
1090  FORMAT(*  FOR  SATURATED  LIQUID  ENTER  "0",  FOR  VAPOR  ENTER  "1"*) 

READ  1010, IP 

PRINT  1210 

IFdV.EQ.DGO   TO   330 

PRINT  1095 
1095  FORMAT(*  TO  ENTER  WITH  TEMPERATURE  ENTER  "0",  FOR  PRESSURE  "1"*) 

READ  1010,11 

PRINT  1210 

IFdI.EQ.DGO  TO  290 
250  IFdU.EO.DGO  TO  260 
1100  FORMAT(*  ENTER  A  TEMPERATURE  IN  DEGREES  F*) 

PRINT  1100 

READ  *,TI 

PRINT  1210 

T=(TI-32.)/l. 8+273. 15 

IF(T.LE.O.O)GO  TO  110 

GO  TO  270 
260  PRINT  1110 
1110  FORMAT (*  ENTER  A  TEMPERATURE(K)*) 

READ  *,T 

PRINT  1210 
270  IF(T.LT..000001)G0  TO  110 

IF(T.GT.TCC.OR.T.LT.TTP)GO  TO  280 

P=VPN(T) 

IF{IP.EO.O)P=P+. 00001 

D=FIND  D(P,T) 

CALL  RE  PRO(P,D,T,IU,IV,IC,IP,TF,DELT) 

GO  TO  250 
280  PRINT  1120, TTP,TCC,(TTP-273.15)*1. 8+32., (TCC-273.15)*1. 8+32. 

<  *,F6.2,*  K*,/, 


1120 

FORMAT(*  FOR  SATURATION  * 

,F6 

.2,*  < 

TEMP 

A  *  OR  *,F7.2,*  < 

TEMP  <  * 

,F7 

.2,*  F*) 

GO  TO  250 

290 

IFdU.EQ.DGO  TO 
PRINT  1130 

300 

1130 

FORMAT(*  ENTER  A 
READ  *,PI 
PRINT  1210 
IF(PI.LE.O.O)GO 
P=(PI/14. 695949)- 
GO  TO  310 

PRESSURE 

TO  110 
M01325 

IN 

LB/SQ 

IN*) 
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300  PRINT  1140 
1140  FORMAK*  ENTER  A  PRESSURE(MPA)*) 
READ  *,PI 
PRINT  1210 

IF(PI.LE.0.0)G0  TO  110 
P=PI 
310  IF{P.GT.PCC.OR.P.LT.PTP)GO  TO  320 
T=FIND  TV(P) 
P=VPN(T) 

IF(IP.E0.0)P=P+. 00001 
D=FIND  D(P,T) 

CALL  RE  PRO(P,D,T,IU,IV,IC,IP,TF,DELT) 
GO  TO  290 
320  PRINT  1150,  PTP,PCC 
1150  F0RMAT(*  YOUR  INPUT  PRESSURE  IS  OUTSIDE  THE  RANGE  OF  SATURATION* 
A  *  PRESSURES*/*  FOR  THIS  FLUID.  TRIPLE  P0INT=*,F6.5 ,*  MPA,* 
B  *  CRITICAL  P0INT=*,F6.3,*  MPA*/*  TRY  AGAIN*) 
GO  TO  290 
330  IFdC.EQ.DGO  TO  370 
IFdU.EO.DGO  TO  340 
PRINT  1160 
1160  F0RMAT(*  ENTER  A  STARTING  TEMPERATURE,  A  FINAL  TEMPERATURE* 

A  /*  AND  A  TEMPERATURE  INCREMENT,  IN  DEGREES  F  AND  IN  THAT  ORDER*) 
READ  *,TS,TF,DELT 
PRINT  1210 

IF(DELT.LE.0.0)G0  TO  110 
TS=(TS-32.)/l. 8+273. 15 

TF=(TF-32.)/l. 8+273. 15  .  • 

DELT=DELT/1.8 

IF(TS.LT.TTP.OR.TS.GT.TCC)GO  TO  360 
IF(TF.LT.TTP.0R.TF.GT.TCC)G0  TO  360 
GO  TO  350 
340  PRINT  1170 
1170  FORMAT!*  ENTER  A  STARTING  TEMPERATURE,  A  FINAL  TEMPERATURE* 
1/*  AND  A  TEMPERATURE  INCREMENT  IN  KELVINS  AND  IN  THAT  ORDER*) 
READ  *,TS,TF,DELT 
PRINT  1210 

IF(DELT.LE.0.0)G0  TO  110 
IF{TS.LT.TTP.OR.TS.GT.TCC)GO  TO  360 
IF(TF.LT.TTP.OR.TF.GT.TCC)Gn  TO  360 
350  T=TS 

P=VPN(T) 

IF(IP.EQ.0.0)P=P+. 00001 
D=FIND  D(P,T) 

CALL  RE  PRO(P,D,T,IU,IV,IC,IP,TF,DELT) 
GO  TO  330 
360  PRINT  1180, TTP,TCC,(TTP-273.15)*1. 8+32.  ,(TCC-273.15)*1. 8+32 
1180  F0RMAT(*  FOR  SATURATION,  *,F6.2,*  <  TEMP  <  *,F6.2,*  K*,/, 
A  ,13X,*0R,  *,F7.1,*  <  TEMP  <  *,F7.1,*  F.  TRY  AGAIN.*) 
GO  TO  330 
370  IFdU.EO.DGO  TO  380 
PRINT  1190 
1190  FORMATf*  ENTER  PRESSURE(PSIA) ,  STARTING  TEMPERATURE(F) ,  FINAL  * 
A  *TEMPERATURE(F)*/*  AND  A  TEMPERATURE  INCREMENT,  IN  THAT  ORDER*) 
READ  *,PI,TS,TF,DELT 
PRINT  1210 

IF(DELT.LE.O.O)GO  TO  110 
P=(PI/14. 695949)*. 101325 
T=(TS-32.)/l. 8+273. 15 
TF=(TF-32.)/l. 8+273. 15 
DELT=DELT/1.8 
CALL  LIMITS(P,T,IL) 
IF(IL.LE.O)GO  TO  370 
CALL  LIMITS(P,TF,IL) 
IF(IL.LE.O)GO  TO  370 
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GO  TO  390 
380  PRINT  1200 
1200  FORMAT(*  ENTER  A  PRESSURE(MPA)  ,  A  STARTING  TEMPERATURE(K)  ,  A 
1*FINAL  TEMPERATURE  AND  A*/*  TEMPERATURE  INCREMENT* 
2*  AND  IN  THAT  ORDER,  TO  RESTART  PROGRAM  ENTER  "0,0,0,0"*) 
READ  *,PI,TS,TF,DELT 
PRINT  1210 

IF(DELT.LE.0.0)G0  TO  110 
T=TS 
P=PI 

CALL  LIMITS(P,T,IL) 
IF(IL.LE.O)GO  TO  370 
CALL  LIMITS(P,TF,IL) 
IF(IL.LE.O)GO  TO  370 
390  D=FIMD  D(P,T) 

CALL  RE  PRO(P,D,T,IU,IV,IC,IP,TF,DELT) 
GO  TO  370 
999  CONTINUE 
1210  FORMAT(*  *) 
END 


SUBROUTINE  REPRO(P,D ,T, lU, IV, IC, IP ,TF ,DELT) 

DIMENSION  G(32)  ,VP(q) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

COMMON/CONT/IF 

COMMON/CRIT/EM,EOK,RM,TC,DC,X,PC,SIG 

C0MM0N/CPID/GI(11),GH(11),GL(11) 

COMMON/HAN/CR,TCI 

N=500 

IF(IV.EO.O)TF=T-1. 

IF(IU.E0.O)GO  TO  100 

PRINT  1000 

PRINT  1010 

PRINT  1020 

GO  TO  110 
100  PRINT  1000 

PRINT  1030 

PRINT  1040 
110  CONTINUE 

DO  210  1=1, N 

IFd.EO.DGO   TO    120 

D^FIND   D(P,T) 
120  H=ENTHAL{P,D,T) 

E=H-1000.*P/D 

S=ENTROP(D,T) 

W=SOUND(D,T) 

CPP=CPfD,T) 

CVV=CV(D,T) 

IF(IU.EQ.O)GO  TO   160 

IF(IF.EQ.4.0R.IF.E0.6)   GO   TO    150 
C  IF-4  FOR   N2,      IF=6  FOR   02. 

IF(IF.En.l)    GO   TO   140 
C  IF=1    FOR   AR 

IF(IF.E0.3)   GO   TO   130 
C  IF=3  FOR   H2 

C  IF=2  FOR   C2H4,    IF-5  FOR   NF3 

PRINT  2000,   P,T,D,E,H,S,CVV,CPP,W 

GO   TO   200 
130   EPS=DIEL(P,D,T) 

PRINT  2010,   P,T,D,E,H,S,CVV,CPP,W,EPS 

GO  TO   200 
140   TH=THERM(D,T)*1000. 

V=VISC(D,T) 

PRINT  2020,    P,T,D,E,H,S,CVV,CPP,W,V,TH 
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GO  TO  200 
150  V=VISC(D,T) 

TH=THERM(D,T)*1000. 
EPS=DIEL(P,D,T) 

PRINT  2030,  P,T,D,E,H,S,CVV,CPP,W,V,TH,EPS 
GO  TO  200 
160  H=H/(2.324445*EM) 
E=E/(2.324445*EM) 
S-S/(4.184001*EM) 
CPP=CPP/{4.184001*EM) 
CVV=CVV/(4.184001*EM) 
W=W*3. 280840 

P0=(P/.101325)*14. 695949 
D0=D*EM/16. 01846371 
T0=T*1. 8-459. 67 

IF{IF.EQ.4.0R.IF.EQ.6)  GO  TO  190 
IF(IF.EO.l)  GO  TO  180 
IF(IF.E0.3)  GO  TO  170 
PRINT  3000,  PO,TO,DO,E,H,S,CVV,CPP,W 
GO  TO  200 
170  EPS=DIEL(P,D,T) 

PRINT  3010,  PO,TO,DO,E,H,S,CVV,CPP,W,EPS 
GO  TO  200 
180  V=VISC(D,T)*. 067196897 
TH=THERM(D,T)*. 578176 

PRINT  3020,  P0,T0,D0,E,H,S,CVV,CPP,W,V,TH 
GO  TO  200 
190  V=VISC(D,T)*. 067196897 
TH-THERM(D,T)*. 578176 
EPS=DIEL(P,D,T) 

PRINT  3030,  P0,TO,D0,E,H,S,CVV,CPP,W,V,TH,EPS 
200  T=T+DELT 

IF(T.GT.TF+.01)G0  TO  220 
IF(IC.EQ.O)P-VPN(T) 
IF(IP.EQ.O)P=P+. 00001 
210  CONTINUE 
220  CONTINUE 
PRINT  1000 
1000  FORMAT(*  *) 

RETURN 
1010  FORMAT (3X,*P*,5X,*T*,6X,*DEN*,5X,*E*,6X,*H*,6X,*S*,4X,*CV*,4X, 

A  *CP* ,3X ,*SOUND* ,2X ,*VISC* ,2X  ,*C0ND*  ,2X  ,*DIEL*) 
1020  F0RMAT(3X,*MPA*,3X,*K*,6X,*M0L/L*,2X,*===  J/HOL  ==*,2X,"====  ", 

A  "J/MOL-K  -==" ,2X ,*M/S* ,3X ,"PA-S" ,1X , "MW/M-K" ,2X ,*==*,/ ,61X ,*E+6*) 
1030  F0RMAT(3X,*P*,5X,*T*,7X,*DENS*,3X,*E*,8X,*H*,5X,*S*,5X,*CV*, 

A  4X,*CP*,2X,*S0UND*,1X,*VISC*,2X,*C0ND*,2X,*DIEL*) 
1040  F0RMAT(3X,*PSIA*,2X,*F*,7X,*LB/*,3X,*-===  BTU/  ===*, IX ,*======*, 

A  *  BTU/  ====-*, 2X,*F/S*,3X,*LB/*,2X,*BTU/*,2X,*=-*,/,17X,*CU  FT*, 

B  6X,*LB*,14X,"LB-F",13X,"FT-S",lX,"FT-HR-F",/,62X,"E+7") 
2000  F0RMAT(F7.3,F8.3,F7.3,2F7.0,F6.1,F5.1,F6.1,F6.0) 
2010  F0RMAT(F7.3,F8.3,F7.3,2F7.0,F6.1,F5.1,F6,1,F6.0,13X,F8.5) 
2020  F0RMAT(F7.3,F8.3,F7.3,2F7.0,F6.1,F5.1,F6.1,F6.0,F6.1,F6.1) 
2030  F0RMAT(F7.3,F8.3,F7.3,2F7.0,F6.1,F5.1,F6.1,F6.0,F6.1,F6.1,F8.5) 
3000  F0RMAT(F7.1,F8.3,F7.3,2F7.1,3F6.3,F6.0) 
3010  F0RMAT(F7.1,F8.3,F7.3,2F7.1,3F6.3,F6.0,11X,F8.5) 
3020  F0RMAT(F7.1,F8.3,F7.3,2F7.1,3F6.3,F6.0,F5.2,F6.4) 
3030  F0RMAT(F7.1,F8.3,F7.3,2F7.1,3F6.3,F6.0,F5.2,F6.4,F8.5) 
END 
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Appendix  E.     Data  Subroutines 


SUBROUTINE   DATA  AR 
DIMENSION   G(32),VP(9) 

DIMENSION   GV(9) ,GT(9) ,FV(4) ,FT(4) ,EV(8) ,ET{8) 
DIMENSION   A(20) 

COMMON/SEN/BETA, XO, DELTA, El,    E2,    AGAM 
COMMON/SATC/A,DTPV 

COMMON/CRIT/  EM,    EOK,    RM,    TC,    DC,  X    ,    PC,    SIG 
COMMON/DATA/G,R , GAMMA, VP , DTP ,PCC ,PTP ,TCC ,TTP ,TUL ,TLL ,PUL ,DCC 
C0MM0N/DATA1/GV,GT,FV,FT,EV,ET 
COMMON/CPID/GKll)  ,GH(11)  ,GL(11) 
COMMON/ISP/N,NW,NWW 
COMMON/FIXPT/TO,SO,HO 
C0MM0N/DIEL/BX(6),PX(6) 
NWW=0 
PRINT   100 
100  FORMAT(*  THE   TEMPERATURE    RANGE   FOR   ARGON    IS  83.8  TO  400   K* 
1/*  WITH   PRESSURES   TO   100  MPA*) 
N=0   $  NW-1 

EM=39.948  $  E0K=152.8  $  RM=3.669E-08  $  TC=150.725 
DC-0.533  $  X=1.7124  $  PC=4.8619  $  SIG=3.297 
X0=0.183     $   BETA=0.355   $  DELTA=4.352  $  El=2.27   $  E2=0.287 


AGAM=1.190 


C 


ARGON  COEFFICIENTS  FOR  MBWR  EON 


G 

[    1) 

G 

2) 

G 

3) 

G 

4) 

G 

5) 

G 

6) 

G 

7) 

G 

8) 

G 

9) 

G 

10) 

G 

11) 

G 

12) 

G 

13) 

G 

14) 

G 

15) 

G 

16) 

G 

17) 

G 

18) 

G 

19) 

G 

20) 

G( 

21) 

G 

22) 

G 

23) 

G 

2^) 

G 

25) 

G 

26) 

G 

27) 

G 

28) 

G 

29) 

G 

30) 

G 

31) 

G 

32) 

GAMMA= 


-.65697312940E-04 

.18229578010E-01 
-.36494701410E+00 

.12320121070E+02 
-.86135782740E+03 

.79785796910E-05 
-.29114891100E-02 

.75818217580E+00 

.87804881690E+03 

.14231459890E-07 

.16741461310E-03 
-.32004479090E-01 

.25617663720E-05 
-.54759349410E-04 
-.45050320580E-01 

.20132546530E-05 
-.16789412730E-07 

.42073292710E-04 
-.54442129960E-06 
-.80048550110E+03 
-.13193042010E+05 
-.49549239300E+01 

.80921321770E+04 
-.98701040510E-02 

.20204415620E+00 
-.16374172050E-04 
-.70389441360E-01 
-.11543245390E-07 

.15559901170E-05 
-.14921785360E-10 
-.10013560710E-08 

.29339632160E-07 
.0055542372 


C  ****  ARGON  COEFFICIENTS  FOR  VAPOR  PRESSURE 
C 

VP(1)=   3.4151115519 


$  VP(2)=   1.1910812519 
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VP(3)=   -.3407632334 
VP(6)=   1.5 

VP(7)=  83.80 
VP(4)=   0.0 


S  VP(5)=  .89555855251 
$  VP(9)=  .06890606625 
$  VP(8)=  150.86 


ARGON  COEFFICIENTS  FOR  SATURATED  LIQUID  AND  VAPOR  DENSITIES, 


1)  = 
Z)  = 

3)  = 

4)  = 

5)  = 
6)- 
7)  = 


A( 

A( 

A( 

A( 

A( 

A( 

A( 

A(  8)  = 

A(  9)  = 

A(10)= 

A(ll)= 

A(12)= 

A(13)= 

A(14)= 

A(15)= 

A(16)= 

A(17)= 

A(18)= 

A(19)= 

A(20)= 


-.270262923777E+02 

.131040241866E+00 
-.267486438128E+01 

.300176804406E+02 
-.875899149326E+02 
-.408267436456E+02 

.104268066451E+03 
-.671278555379E+02 

.151002935701E+02 
-.331243536637E+02 

.633146212581E+02 
-.427149706899E+02 

.100599900030E+02 

.137682084900E+02 
-.664630363191E-01 

.133368782730E+01 
-.144371463244E+02 

.501938472000E+02 
-.230888463887E+02 

.465318358887E+01 


****  ARGON  COEFFICIENTS  FOR  IDEAL  GAS  CP. 

GI(1)=GI(2)-GI(3)=0 

GI(4)=GI(9)-2.5 

GI(5)=GI(6)=GI(7)=0 

GI(8)=0 

GI(10)=GI(11)=0.0 

****  ARGON  COEFFICIENTS  FOR  FIRST  TERM  OF  VISCOSITY  EQN. 
****  VISCOSITY  IN  MICRO-PA*S,  DENSITY  IN  G/CC. 

GV(1)=  .61145472787E+04  S  GV(2)=  -.10394390312E+05 

GV(3)=   .67594614619E+04  S  GV(4)=  -  .22536509380E+04 

6V(5)-  .42593950138E+03  S  GV(6)-  -.47252671093E+02 

GV(7)=   .31795275425E+01  S  GV(8)=  -  .11629083780E+00 
GV(9)-  .18043010592E-02 

****  ARGON  COEFFICIENTS  FOR  FIRST  TERM  OF  THERMAL  CONDUCTIVITY  EQN. 
****  THERMAL  CONDUCTIVITY  IN  WATT/(M*K),  DENSITY  IN  G/CC. 

GT(1)-  .62777703742E+01  S  GT(2)=  -.96096376637E+01 

GT(3)-  .58887549191E+01  S  GT(4)=  -  .18920926320E+01 

GT(5)=  .34886571437E+00  S  GT(6)=  -.38016786193E-01 

GT(7)=  .25207283167E-02  $  GT(8)=  -.91098744478E-04 
GT(9)=  .13990842942E-05 

****  ARGON  COEFFICIENTS  FOR  SECOND  TERM  OF  VISCOSITY  EQN. 

FV(1)=   .14653652433E+00  S  FV(2)=  -.77487424965E-01 

FV(3)=   .14000000000E+01  S  FV(4)=  .15280000000E+03 

****  ARGON  COEFFICIENTS  FOR  SECOND  TERM  OF  THERMAL  CONDUCTIVITY  EQN, 

FT(1)-   .24142103270E-01  S  FT(2)=  .75696234255E-02 

FT(3)=   .lOOOOOOOOOOE+01  $  FT(4)-  .15280000000E+03 

****  ARGON  COEFFICIENTS  FOR  THIRD  TERM  OF  VISCOSITY  EQN. 
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EV(1)=  -.12313579086E+02  $  EV(2)=  .20694685712E+00 
EV(3)=  .16029145122E+02  $  EV(4)=  .11717461351E+04 
EV(5)=  -.56995898780E+03  $  EV(6)=  .40136071933E+02 
EV(7)=  .39870122403E+02  $  EV(8)=   .53700000000E+00 

C 

C  ****  ARGON  COEFFICIENTS  FOR  THIRD  TERM  OF  THERMAL  CONDUCTIVITY  EQN. 

C 

ET(1)=  -.33327027332E+02  $  ET(2)=  0.0 

ET(3)=   .30694859971E+02  $  ET(4)=  0.0 

ET(5)=  .22956551674E+04  $  ET(6)=  -.35559415848E+03 

ET(7)=  0.0  $  ET(8)=  1.0 

C 

C  ****  ARGON  COEFFICIENTS  FOR  MELTING  EQUATION. 

C 

PX(1)=  -.210562165E+03 
PX(2)=  .177760527E+00 
PX(3)=  .159817868E+01 

C     THE  FOLLOWING  ARE  CONSTANTS  USED  IN  THE  COMPUTATION  OF  PROPERTIES  FOR 
C     ARGON 

R=. 08205616*. 101325 

T0=298.15  $  S0=154.7335  $  H0=6196.5 

GI(10)-H0-HI(T0) 

GI(11)=S0-SI(T0) 

DTP=.3540027619188E+02 

DTPV=. 1029227022329 

TCC=VP(8) 

PCC=VPN(TCC) 

PTP=VP(9) 

TTP=VP(7) 

TUL=400. 

TLL=TTP 

PUL=100. 

DCC=13.41 

RETURN 

END 


C 
C 

SUBROUTINE  DATA  C3 
C     COEFFICIENTS  FOR  ETHYLENE 

DIMENSION  A(20) 

DIMENSION  G(32),VP(9) 

DIMENSION  GV(9),GT(9),FV(4),FT(4),EV(8),ET(4) 

C0MM0N/CPID/GI(11),GH(11),GL(11) 

COMMON/CRIT/  EM,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

C0MM0N/DATA1/GV,GT,FV,FT,EV,ET 

COMMON/ISP/N,NW,NWW 

COMMON/SATC/A,DTPV 

C0MM0N/DIEL/BX(5),PX(6) 

N=0 

PRINT  100 
100  FORMAT(*  THE  TEMPERATURE  RANGE  FOR  ETHYLENE  IS  104  TO  400K* 
1/*  WITH  PRESSURES  TO  40  MPA*) 
C 

C  ****  ETHYLENE  COEFFICIENTS  FOR  MBWR  EQN 
C 

G(  1)=  -.2146684366683E-02 

G(  2)=  .1791433722534E+00 

G(  3)=  -.3675315603930E+01 

G(  4)=  .3707178934669E+03 

G(  5)=  -.3198282566709E+05 

G(  6)=  .5809379774732E-04 

G(  7)=  -.7895570824899E-01 
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8)=  .1148620375835E+02 
9)=  .2713774629iq3E+05 
10)=  -.8647124319107E-05 
11)=  .1617727266385E-01 
12)=  -.2731527496271E+01 
13)=  -.2672283641459E-03 
14)=  -.4752381331990E-02 
15)=  -.6255637346217E+01 
16)=  .4576234964434E-03 
17)=  -.7534839269320E-05 
18)=  .1638171982209E-01 
19)=  -.3563090740740E-03 
20)=  -.1833000783170E+05 
21)=  -.1805074209985E+07 
22)=  -.4794587918874E+03 
23)=  .3531948274957E+07 
24)=  -.2562571039155E+01 
25)=  .1044308253292E+03 
26)=  -.1695303363659E-01 
27)=  -.1710334224958E+03 
28)=  -.2054114462372E-04 
29)=  .6727558766661E-02 
30)=  -.1557168403328E-06 
31)=  -.1229814736077E-04 
32)=  .4234325938573E-04 
GAMMA=-.0172 


****  ETHYLENE  COEFFICIENTS  FOR  VAPOR  PRESSURE 


VP(1)=   8.2095798 
VP(3)=  -1.692585975 
VP(6)=   1.5 
VP(7)=  103.986 
VP(4)=   -.1976495575 


VP(2)= 
VP(5)= 
VP(9)= 


4.315424145 
3.446501098 
.00012129514 


$  VP(8)=  282.3428 


ETHYLENE  COEFFICIENTS  FOR  SATURATED  LIQUID  AND  VAPOR  DENSITIES. 

1)=  -.609621515594E+02 

2)=  .203185312702E-01 

3)=  -.925441265813E+00 

4)=  .243630795888E+02 

5)=  -.854745622888E+03 

6)=  .123927868183E+04 

7)=  -.142710711789E+04 

8)=  .837358670405E+03 

9)=  .432203696552E+03 

10)=  -.137917541161E+04 

11)=  .126858600124E+04 

12)=  -.571552321713E+03 

13)=  .106012234360E+03 

14)=  -.479047060183E+01 

15)=  .151381345283E-01 

16)=  -.403456079445E+00 

17)=  .508683920225E+01 

18)=  -.246711997987E+02 

19)=  .980030915247E+01 

20)=  -.216846516122E+01 


****  ETHYLENE  COEFFICIENTS  FOR  IDEAL  GAS  HEAT  CP , 

GI(l)=.5603615762E+6 

GI(2)=-.2141069802E+5 

GI(3)=.2532008897E+3 

GI(4)=.3554495281E+1 

GI(5)=-.9951927478E-2 
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GI(6)=.5108931070E-4 

GI(7)=-.1928667482E-7 

GI(8)=-.2061703241E+2 

GI(9)=.3E+4 

GI(10)=GI(11)=0 

T0=298.15  $  50=219.451   $  H0=10564.6 

GI(10)=H0-HI(T0) 

GI(11)=S0-SI(T0) 
C 

C  ****  ETHYLENE  COEFFICIENTS  FOR  MELTING  EQUATION. 
C 

PX(1)=  -.357923875E+03 

PX(2)=  .245332143E-01 

PX(3)=  .206450000E+01 

C     THE  FOLLOWING  CONSTANTS  USED  IN  THE  COMPUTATION  OF  PROPERTIES  FOR 
C     ETHYLENE 

GAMMA=-.0172 

R=. 08205616*. 101325 

DTP=23. 34296694034 

DTPV=.1425455127094E-03 

EM=28.054 

TCC=VP(8) 

PCC=VPN(TCC) 

PTP=VP(9) 

TTP=VP(7) 

TUL=400. 

TLL=TTP 

PUL=40. 

DCC=7.650 

RETURN 

END 


SUBROUTINE  DATA  P  H2 

■DIMENSION  G(32),VP(9),GI(11),GH(11),GL(11) 

DIMENSION  GV(9),GT(9),FV(4),FT(4),EV(8),ET(8) 

DIMENSION  A(20) 

COMMON/SATC/A,DTPV 

COMMON/CPID/GI,GH,GL 

COMMON/CRIT/  EM,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

C0MM0N/DATA1/GV,GT,FV,FT,EV,ET 

C0MM0N/DIEL/BX(6),PX(6) 

COMMON/ISP/N,NW,NWW 

N=l 

PRINT  100 
100  FORMAT(*  THE  TEMPERATURE  RANGE  FOR  HYDROGEN  IS  13.8  TO  400  K* 
1/*  WITH  PRESSURES  TO  120  MPA*) 
C 

C  ****  PARAHYDROGEN  COEFFICIENTS  FOR  MBWR  EQUATION. 
C 

G(  1)=  .4675528393416E-04 

G(  2)=  .4289274251454E-02 

G(  3)=  -.5164085596504E-01 

G(  4)=  .2961790279801E+00 

G(  5)=  -.3027194968412E+01 

G(  6)=  .1908100320379E-05 

G(  7)=  -.1339776859288E-03 

G(  8)=  .3056473115421E-01 

G(  9)-  .5161197159532E+01 

G(10)=  .1999981550224E-07 

G(ll)=  .2896367059356E-04 

G(12)=  -.2257803939041E-02 

G(13)=  -.2287392761826E-06 

G(14)=  .2446261478645E-05 
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G(15)=  -.1718181601119E-03 

G(16)=  -.5465142603459E-07 

G(17)=  .4051941401315E-09 

G(18)=  .1157595123961E-06 

G(19)=  -.1269152728389E-08 

G(20)=  -.4983023605519E+01 

6(21)=  -.1606676092098E+02 

G(22)=  -.1926799185310E-01 

G(23)=  .9319894638928E+00 

G(24)=  -.3222596554434E-04 

G(25)=  .1206839307669E-03 

G(26)=  -.3841588197470E-07 

G(27)=  -.4036157453608E-05 

G(28)=  -.1250868123513E-10 

G(29)=  .1976107321888E-09 

G(30)=  -.2411883474011E-13 

G(31)=  -.4127551498251E-13 

G(32)=  .8917972883610E-12 
GAMMA=-.0041 

'  PARAHYDROGEN  COEFFICIENTS  FOR  VAPOR  PRESSURE 


VP(1)=   3.05300134164  $  VP(2)= 

VP(3)=  -0.655461216567  $  VP(5)= 

VP(6)=   1.5814454428  $  VP(9)= 

VP(7)=  13.8  $  VP(8)= 
VP(4)=   0.0 


2.80810925813 
1.59514439374 
0.0070420875 
32.938 


PARAHYDROGEN  COEFFICIENTS.   SATURATED  LIQUID  AND  VAPOR  DENSITIES, 

1)=   .916617720187E+02 

2)=  -.1794Q2524446E+00 

3)=   .454671158395E+01 

4)=  -.658499589788E+02 

5)=   .734466804535E+03 

-6)=  -.682501045175E+03 

7)=   .631783674710E+03 

8)=  -.539408873282E+03 

9)=   .430923811783E+03 

10)=  -.300295738811E+03 
11)=   .156567165346E+03 

12)=  -.504103608225E+02 
13)=   .720706926514E+01 

14)=  -.123944440318E+03 
15)=   .140334800142E+01 

16)=  -.211023804313E+02 
17)=   .173254622817E+03 

18)=  -.444294580871E+03 
19)=   .138699365355E+03 

20)=  -.235774161015E+02 


****  PARAHYDROGEN  COEFFICIENTS  FOR  IDEAL  GAS  CP. 
FOR  TEMPERATURES  >  140  K. 


GH(  1)= 

.5262185164597E+08 

GH(  2)= 

-.1487906248823E+07 

GH(  3)= 

.1601391392264E+05 

GH(  4)= 

-.8031235938946E+02 

GH(  5)= 

.2307407941873E+00 

GH(  6)= 

-.3176386248370E-03 

GH(  7)= 

.1643857271214E-06 

GH(  8)= 

.9230816464058E+01 

GH(  9)= 

.3E+4 

GH(10)= 

0.0 

GH(11)= 

0.0 
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c 

C  ****  PARAHYDROGEN  COEFFICIENTS  FOR  IDEAL  GAS  CP. 

C      FOR  TEMPERATURES  40  <  T  <  140  K. 

C 

GL(  1)=  .2905965792270E+06 

GL(  2)=  -.2831103639248E+05 

GL(  3)=  .1050424877391E+04 

GL(  4)=  -.1535751501769E+02 

GL(  5)=   .1218941696566E+00 

GL(  6)=  -.2599406479908E-04 

GL(  7)=  -.1288757333406E-05 

GL(  8)=  .1717441975231E+06 

GL(  9)=  .3E+4 

GL(10)=  0.0 

GL(11)=  0.0 

T0=298.15  $  S0=130.407   $  H0=8409.8 
C        HO  AND  SO  TAKEN  FROM  WOOLLEY,  SCOTT,  AND  BRICKWEDDE 
C        NBS  RP  1932,  PAGE  387.  VALUES  FOR  PARAHYDROGEN. 

GH(10)=H0-HI(T0) 

GL(10)=GH(10) 

GH(11)=S0-SI(T0) 

GL(11)=GH(11) 
C 
C 

C  ****  PARAHYDROGEN  COEFFICIENTS  FOR  MELTING  PRESSURE. 
C      FOR  TEMPERATURES  <  22  K  : 
C 

PX(1)=  -21.281484395 

PX(2)-    .125746643 

PX(3)=   1.955 
C 

C  ****  PARAHYDROGEN  COEFFICIENTS  FOR  MELTING  PRESSURE. 
C      FOR  TEMPERATURES  >  22  K  : 
C 

PX(4)=  -26.5289115 

PX(5)=    .248578596 

PX(6)=   1.764739 
C 

C  ****  PARAHYDROGEN  COEFFICIENTS  FOR  DIELECTRIC  CONSTANT. 
C 

BX(1)=  .20245443E-02 

BX(2)=  .37171832E-06 

BX(3)=  -.92085013E-08 

BX(4)=  -.34065328E-11 

BX(5)=  .0 

BX(6)=   .0 

C     THE  FOLLOWING  ARE  CONSTANTS  USED  IN  THE  COMPUTATION  OF  PROPERTIES  FOR 
C     PARAHYDROGEN 

R=. 08205616*. 101325 

DTP=.3821428945438E+02 

DTPV=.6322296353698E-01 

EM=2. 01594 

TCC=VP(8) 

PCC=VPN(TCC) 

PTP=VP(9) 

TTP=VP(7) 

TUL=400. 

TLL=TTP 

PUL=120. 

DCC=15.556 

RETURN 

END 
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SUBROUTINE  DATA  N2 

DIMENSION  G(32),VP(9) 

DIMENSION  GV(9),GT(9),FV(4),FT(4),EV(8),ET(8) 

DIMENSION  A(20) 

COMMON/CRIT/  EM,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 

COMMON/SATC/A,DTPV 

C0MM0N/DATA1/GV,GT,FV,FT,EV,ET 

COMMON/SEN/BETA, XO, DELTA, El,  E2,  AGAM 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

C0MM0N/CPID/GI(11),GH(11),GL(11) 

COMMON/ISP/N,NW,NWW 

C0MM0N/DIEL/BX(6),PX(6) 

NWW=0 

PRINT  100 
100  FORMAT(*  THE  RANGE  OF  TEMPERATURE  FOR  NITROGEN  IS  63.15  TO  1900K* 
1/*  WITH  PRESSURES  TO  1000  MPA*) 

N=0  $  NW=1 

E0K=118.  $  RM=3.933E-08  $  TC=126.24 
X=l. 67108  l;  PC=3.443  $  SIG=3.54 
BETA=0.355  $  DELTA=4.352  $  El=2.17  S  E2=0.287 


EM=28.016  : 
DC=0.3139  3 
X0=0.164  ; 
AGAM=1.190 


NITROGEN  COEFF  FOR  MBWR  EQUATION 

1)=  .138029747465691E-03 
2)=  .108450650134880E-01 
3)=  -.247132406436209E+00 
4)=  .345525798080709E+01 
5)=  -.427970769066595E+03 
6)=   .106491156699760E-04 
7)=  -.114086707973499E-02 
8)=  .144490249728747E-04 
9)=  .187145756755327E+04 
10)=  .821887688683079E-08 
11)=  .236099049334759E-03 
12)=  -.514480308120135E-01 
13)=  .491454501366803E-05 
14)=  -.115162716239893E-03 
15)=  -.716803724664983E-01 
16)=   .761666761949981E-05 
17)=  -.113093006621295E-06 
18)=  .373683116683089E-04 
19)=  -.203985150758086E-06 
20)=  -.171966200898966E+04 
21)=  -.121305519974777E+05 
22)=  -.988139914142789E+01 
23)=  .561988689351085E+04 
24)=  -.182304396411845E-01 
25)=  -.259982649847705E+00 
26)=  -.419189342315742E-04 
27)=  -.259640667053023E-01 
28)=  -.125868320192119E-07 
29)=   .104928659940046E-05 
30)=  -.545836930515201E-10 
31)=  -.767451167059717E-09 
32)=  .593123287099439E-08 
GAMMA=-.0056 


****  NITROGEN  COEFF  FOR  VAPOR  PRESS. 


VP(1)=  5.1113192094 

VP(3)=  -.15108730916 

VP(6)=  1.5 

VP(7)=  63.15 


VP(2)=. 6482667539 
VP(5)=. 74028493342 
VP(9)=. 012462975 
VP(8)=  126.26 
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VP(4)=   0.0 
'     HZ   COEF  FOR  SAT  LIQUID  AND  VAPOR  DENS. 


A(  1)=  -.158453465507E+02 

A(  2)=  .419136911423E-01 

A(  3)=  -.101965371660E+01 

A(  4)=  .134763743799E+02 

A(  5)=  -.109930399087E+03 

A(  6)=  .925518835497E+02 

A(  7)=  -.956233831320E+02 

A(  8)=  .100104366710E+03 

A(  9)=  -.701857937398E+02 

A(10)=  .900076998647E+01 

A(ll)=  .286981120347E+02 

A(12)=  -.216767601780E+02 

A(13)=  .496558226471E+01 

A(14)=  .218307928477E+02 

A(15)=  -.126493309807E+00 

A(16)=  .241544188633E+01 

A(17)=  -.245256871794E+02 

A(18)=  .935925207124E+02 

A(19)=  -.36Q938251632E+02 

A(20)=  .757453271989E+01 

****  N2  COEF  FOR  IDEAL  GAS  CP. 

GI(1)=  -0.735210401157252E  03 

GI(2)=  0.342239980411978E  02 

6I(3)=  -0.557648284567620E  00 

GI(4)=  0.350404228308756E  01 

GI(5)=  -0.173390185081005E-04 

GI(6)=  0.174650849766463E-07 

GI(7)=  -0.356892033544348E-11 

GI(8)=  0.100538722808834E  01 

GI(9)=  0.335340610000000E  04 
GI(10)=GI(11)=0.0 

-***  N2  COEF  FOR  1ST  TERM  OF  VISC 

Vise  IN  MICRO  PA*S,  DENSITY  IN  G/CC 

GV(1)=  -.18224240000E+05  $  GV(2)= 

GV(3)=  -.91542324494E+04  $  GV(4)= 

GV(5)=  -.36307214228E+03  $  GV(6)= 

GV(7)=  -.22261880817E+01  $  GV(8)= 
GV(9)=  -.11894029104E-02 


.19915327374E+05 
.23255484059E+04 
.36457506811E+02 
.78053904895E-01 


****  N2  COEF  FOR  1ST  TERM  THERM.  COND. 

THERM.  COND.  IN  WATT/(M*K),  DENS.  IN  G/CC, 


GT(1)=  -.20029573972E+02  $ 

GT(3)=  .80188959378E+01  $ 

GT(5)=  .15363738965E+01  $ 

GT(7)=  .19360547346E-01  $ 
GT(9)=  .15564670935E-04 

****  N2  COEF  FOR  2ND  TERM  VISC. 


GT(2)= 
GT(4)= 
GT(6)= 
GT(8)= 


FV(1)=  -.11217739623E+00  $  FV(2)= 
FV(3)=   .14000000000E+01  $  FV(4)= 

****  N2  COEF  FOR  2ND  TERM  THERM.  COND. 

FT(1)=  .53875666637E-01  $  FT(2)= 
FT(3)=  .12000000000E+01  $  FT(4)= 


.49765746684E+01 
-.55022716888E+01 
-.22974737257E+00 
-.85677385768E-03 


.32912317244E-01 
.11800000000E+03 


.61027911104E-02 
.11800000000E+03 
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****  N2  COEF  FOR  3RD.  TERM  VISC. 

EV(1)=  -.12128154129E+02  $  EV{2)= 

EV(3)=   .16094611148E+02  $  EV(4)= 

EV(5)=  -.80889801180E+03  $  EV(6)= 

EV(7)=  -.21241135912E+01  $  EV(8)= 


.57156092139E+00 
.36954086158E+04 
.68464435640E+02 
.31500000000E+00 


N2  COEF  FOR  3RD  TERM  OF  THERM.  COND. 


ET(1)=  -.38613291627E+02  : 

ET(3)=  .37201743333E+02  : 

ET(5)=  -.39013509079E+02  : 

ET(7)=  .0  ; 

****  N2  COEF  FOR  MELTING. 

PX(1)=  -.160000281E+03 
PX(2)=  .938575502E-01 
PX(3)=  .179500000E+01 

****  N2  COEF  FOR  DIEL.  CONST. 


ET(2)=  .0 

ET(4)=   .0 

ET(6)=  -.31826109485E+02 

ET(8)=  1.0 


BX(1)^ 
BX(2)^ 
BX(3)-  0. 
BX(4)= 
BX(5)  = 
BX(6)  = 


.43993836E-02 
.18932096E-05 

.31450178E-08 
.28592703E-06 
.44566034E-07 


THE  FOLLOWING  ARE  CONSTANTS  USED  IN  THE  COMPUTATION  OF  PROPERTIES  FOR 
NITROGEN 


R=(8.20539E-2)*. 101325 

T0=298.15  $  S0=191.502 

GI(10)=H0-HI(T0) 

GI(11)=S0-SI(T0) 

DTP=.3097717741477E+2 

DTPV=.242822085710E-1 

TCC-VP(8) 

PCC=VPN(TCC) 

PTP=VP(9) 

TTP=VP(7) 

TUL=1900. 

TLL=TTP 

PUL=1000. 

DCC=11.21 

RETURN 

END 


$  H0=8669.0 


SUBROUTINE  DATA  NF3 

DIMENSION  G(32),VP(9),A(20) 

COMMON/CRIT/  EM,  EOK,  RM,  TC ,  DC,  X  ,  PC,  SIG 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

COMMON/SATC/A,DTPV 

C0MM0N/CPID/GI(11),GH(11),GL(11) 

C0MM0N/DIEL/BX(6),PX(6) 

PRINT  100 
100  F0RMAT(*  THE  RANGE  OF  TEMPERATURE  FOR  NITROGEN  TRIFLUORIDE  IS* 
A  ,/,*  66.36  K  TO  500  K,  WITH  PRESSURES  TO  50  MPA.*) 
C 

C  ****  NITROGEN  TRIFLUORIDE  COEFFICIENTS  FOR  MBWR  EQUATION. 
C 

G(  1)=  .1774353868E-02 

G(  2)=  -.5409379418E-01 

G(  3)=   .3976634466E+00 

G(  4)=  -.5209476694E+02 

G(  5)=  -.3285322888E+04 
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6)=  -.5990517411E-04 

7)=  .9217525601E-01 

8)=  -.4848977075E+02 

9)=  -.4235892691E+06 

10)=  -.9824248063E-06 

11)=  .5432235989E-02 

12)=  -.1462388500E+01 

13)=  -.3366180440E-03 

14)=  .2801374599E-01 

15)=  .8435288597E+00 

16)=  -.1324421452E-02 

17)=  .1875604377E-04 

18)=  .2959643991E-01 

19)=  -.7009976870E-03 

20)=  .4365820912E+06 

21)=  -.1111397536E+07 

22)=  .2411866612E+04 

23)=  .3179136276E+06 

24)=  .6166849090E+01 

25)=  .4260854720E+01 

26)=  .1090598789E-01 

27)=  -.3340951059E+01 

28)=  .8597429644E-05 

29)=  .1240544214E-03 

30)=  .1286224248E-07 

31)=  -.8941104276E-07 

32)=  .3353054595E-05 
GAMMA=-.0056 

'  NITROGEN  TRIFLUORIDE  COEFFICIENTS  FOR  VAPOR  PRESSURE, 


VP(1)= 
VP(3)= 
VP(6)= 
VP(7)= 
VP(4)= 


11.593879492 
-2.8727732815 

1.5 
66.36 
-1.37977007 


$  VP(2)=   9.6548502312 

$  VP(5)=   7.3112441673 

$  VP(9)=    .18537827E-6 

$  VP(8)=  234.0 


NITROGEN  TRIFLUORIDE  COEFF.  FOR  SATURATED  LIQUID  AND  VAPOR  DENSITIES 

1)  =  .131285181636E+03 

2)  =  -.226998555536E+01 

3)  =  .384708389498E+02 

4)  =  -.309827268239E+03 

5)  =  -.635348526635E+03 

6)  =  .926979028357E+03 

7)  =  .371498011259E+03 

8)  =  -.182291654470E+04 

9)  =  .105592403853E+04 

10)  =  .121673895344E+04 

11)  =  -.207208928323E+04 

12)  =  .113463394710E+04 

13)  =  -.226642137140E+03 

14)  =  .105087295173E+01 

15)  =  -.131686474246E-02 

16)  =  .390141331900E-01 

17)  =  -.587796597975E+00 

18)  =  .400900521017E+01 

19)  =  -.203385386977E+01 

20)  =  .400434364424E+00 


****     NITROGEN   TRIFLUORIDE   COEFF   CIENTS   FOR    IDEAL  GAS   CP , 

GI(1)=  .7427518245951E+06 
GI(2)=  -.4389825372134E+05 
GI(3)=        .1012629224351E+04' 
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GI(4)=  -.7140693612211E+01 

GI(5)=   .5481339146452E-01 

GI(6)=  -.7677196006769E-04 

GI(7)=   .4203630864340E-07 

GI(8)=  -.6328752997967E+00 

GI(9)=   .3E+4 

GI(10)=GI(11)=0.0 
C 

C  ****  NITROGEN  TRIFLUORIDE  COEFFICIENTS  FOR  MELTING. 
C 

PX(1)=  -.190939971E+03 

PX(2)=  .813750194E-01 

PX(3)=  .185000000E+01 

C     THE  FOLLOWING  ARE  CONSTANTS  USED  IN  THE  COMPUTATION  OF  PROPERTIES  FOR 
C     NITROGEN  TRIFLUORIDE 

C**-   THE  TO, HO, SO  IDEAL  GAS  REFERENCE  VALUES  ARE  FROM 

C  TO, HO, SO  ARE  FROM  NBS  TN  270-3('68),  "SELECTED  VALUES  OF  CHEMICAL 

C  THERMODYNAMIC  PROPERTIES"  WAGMAN  ET  AL.,  PP  70 

T0=298.15  $  S0=260.621   $  HO-11828. 

GI(10)=H0-HI(T0) 

GI(11)=S0-SI(T0) 

EM=71.019 

R=. 0820568*. 101325 

PTP-.185425E-06 

DTP=26.32 

DTPV-.33612E-08 

TTP=66.36 

PCC=4.4607 

DCC=7.92 

TCC=234.0 

TC=234.0 

VP(8)=TC 

DC=7.92 

PUL=  50. 

TUL=500. 

TLL=TTP 

RETURN 

END 


SUBROUTINE  DATA  02 

DIMENSION  G(32),VP(9),GI(11),GH(11),GL(11) 

DIMENSION  GV(9),GT(9),FV(4),FT(4),EV(8),ET(8) 

DIMENSION  A(20) 

COMMON/SATC/A,DTPV 

COMMON/DATAl /GV ,GT , FV , FT , EV , ET 

COMMON/SEN/BETA, XO, DELTA, El,  E2,  AGAM 

COMMON/CRIT/  EM,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

COMMON/ISP/N,NW,NWW 

COMMON/CPID/GI,GH,GL 

C0MM0N/DIEL/BX(6),PX(6) 

N=0  $  NW=1  $  NWW=0 

PRINT  100 
100  FORMAT!*  THE  TEMPERATURE  RANGE  FOR  OXYGEN  IS  54.359  TO  400  K* 
1/*  WITH  PRESSURES  TO  120  MPA*) 

X0=0.183  $  BETA=0.355  $  DELTA=4.352  $  El=2.21  $  E2=0.287 
EM=31.9988  $  E0K=113.0  $  RM=3.8896E-08$  TC=154.575 

DC=0.4362  $  X=2. 210636  $  PC=5.0429   $  SIG-3.437 

AGAM  =  1.190 
C 

C  ****  OXYGEN  COEFFICIENTS  FOR  MBWR  EQUATION. 
C 

G(  1)=  -.4365859650E-04 

G(  2)=   .2005820677E-01 

G(  3)=  -.4197909916E+00 

G(  4)=   .1878215317E+02 
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G(  5)=  -.1287473398E+04 

G(  6)=  .1556745888E-05 

G(  7)=  .1343639359E-03 

G(  8)=  -.2228415518E+00 

G(  9)=  .4767792275E+03 

G(10)=  .4790846641E-07 

G(ll)=  .2462611107E-03 

G(12)=  -.1921891680E-01 

G(13)=  -.6978320847E-06 

G(14)=  -.6214145909E-04 

G(15)=  -.1860852567E-01 

G(16)=  .2609791417E-05 

G(17)=  -.2447611408E-07 

G(18)=  .1457743352E-04 

G{19)=  -.1726492873E-06 

G(20)=  -.2384892520E+03 

G(21)=  -.2301807796E+05 

G(22)-  -.2790303526E+01 

G(23)=  .9400577575E+04 

G(24)=  -.4169449637E-02 

G(25)=  .20n8497853E+00 

G(26)=  -.1256076520E-04 

G(27)=  -.6406362964E-01 

G(28)=  -.2475580168E-08 

G(29)=  .1346309703E-05 

G(30)=  -.116150247nE-10 

G(31)=  -.1034699798E-08 

G(32)-  .2365936964E-07 
GAMMA=-.0056 

'     OXYGEN  COEFFICIENTS  FOR  VAPOR  PRESSURE, 


VP(1)= 
VP(3)= 
VP(6)= 
VP(7)- 
VP(4)= 


7.568956 
-2.13746 

1.514 
54.359 

0.0 


VP(2)=   5.004836 
VP(5)=   3.454481 
VP(9)=    .0001479953 
VP(8)=  154.581 


OXYGEN  COEFFICIENTS  FOR  SATURATED  LIOUID  AND  VAPOR  DENSITIES. 

1)=  .581394753076E+02 

2)=  -.490241196133E-01 

3)=  .168328893252E+01 

4)=  -.325161223398E+02 

5)=  .550300989872E+03 

6)=  -.510968506115E+03 

7)=  .315091559049E+03 

8)=  -.232566659258E+02 

9)=  -.488425479359E+02 

10)=  -.150624217523E+03 

11)=  .280441603851E+03 

12)=  -.176693896861E+03 

13)=  .403247747449E+02 

14)=  .252198688365E+01 

15)=  -.136098316472E-01 

16)=  .282316159403E+00 

17)=  -.286645905341E+01 

18)=  .617024212284E+01 

19)=  -.810220795462E+00 

20)=  -.279601068969E+00 


****  OXYGEN  COEFFICIENTS  FOR  IDEAL  GAS  CP. 

GI(1)=   -0.498199853711943E  04 
GI(2)=    0.230247779995218E  03 
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GI(3)=  -0.345565323510732E  01 

GI(4)=  0.352187677367116E  01 

GI(5)=  -0.435420216024420E-04 

GI(6)=  0.134635345013162E-07 

GI(7)=  0.162059825959105E-10 

GI(8)=  0.103146851572565E  01 

GI(9)=  0.223918105000000E  04 
GI(10)=GI(11)=0.0 

****  OXYGEN  COEFFICIENTS  FOR  FIRST  TERM  OF  VISCOSITY  EQN, 
VISCOSITY  IN  MICRO-PA*S,  DENSITY  IN  G/CC. 


GV(1)=  -.97076378593E+04 
GV(3)=  -.24668758803E+04 
GV(5)=  .37851049522E+02 
GV(7)=  .11134441304E+01 
GV(9)=  .10279379641E-02 


$  GV(2)=  .82801254201E+04 

$  GV(4)=  .21324360243E+03 

$  GV(6)=  -.10487216090E+02 

$  GV(8)=  -.53676093757E-01 


****  OXYGEN  COEFFICIENTS  FOR  FIRST  TERM  OF  THERMAL  CONDUCTIVITY  EQN, 

THERMAL  CONDUCTIVITY  IN  WATTS/(M*K),  DENSITY  IN  G/CC. 

GT(1)=  -.20395052193E+03  $  GT(2)=  .24088141709E+03 

GT(3)=  -.12014175183E+03  $  GT(4)=   .32954949190E+02 

GT(5)=  -.54244239598E+01  $  GT(6)=  .54734865540E+00 

GT(7)=  -.32854821539E-01  $  GT(8)=   .10753572103E-02 
GT(9)-  -.14610986820E-04 

****  OXYGEN  COEFFICIENTS  FOR  SECOND  TERM  OF  VISCOSITY  EQN. 

FV(1)=  .43526515153E+00  $  FV(2)=  -.20361263878E+00 

FV(3)=   .14000000000E+01  $  FV(4)=  .lOOOOOOOOOOE+03 

****  OXYGEN  COEFFICIENTS  FOR  SECOND  TERM  OF  THERMAL  CONDUCITIVY  EQN, 

FT(1)=  .30600000000E-01  $  FT(2)=  .27850000000E-01 

FT(3)-  .11200000000E+01  $  FT(4)=   .lOOOOOOOOOOE+03 

****  OXYGEN  COEFFICIENTS  FOR  THIRD  TERM  OF  VISCOSITY  EQN. 

EV(1)=  -.14454972110E+02  $  EV(2)=  -.31421728994E+00 

EV(3)=   .18201161468E+02  $  EV(4)=   .27390429525E+03 

EV(5)=  -.27498956948E+04  $  EV(6)=  .24340689667E+03 

EV(7)=   .11911504104E+03  $  EV(8)=   .43500000000E+00 

****     OXYGEN  COEFFICIENTS  FOR  THIRD  TERM  OF  THERMAL  CONDUCTIVITY  EQN. 


.0 
.0 
.47350508788E+03 


ET(1)=  -.21520741137E+02  $  ET(2)= 

ET(3)=   .16799504261E+02  $  ET(4)= 

ET(5)=  -.29944878721E+04  $  ET(6)= 

ET(7)=  .0          $  ET(8)=  1.0 

'     OXYGEN  COEFFICIENTS  FOR  MELTING 

PX(1)=  -.267226854E+03 
PX(2)=  .227606348E+00 
PX(3)=  .176900000E+01 


OXYGEN  COEFFICIENTS  FOR  DIELECTRIC  CONSTANT. 


BX(1)= 
BX(2)= 
BX(3)= 


.39608100E-02 
.29700000E-06 
.41300000E-07 


BX(4)=  0.0 
BX(5)=  0.0 
BX(6)=  -.21400000E-07 
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C     THE  FOLLOWING  ARE  CONSTANTS  USED  IN  THE  COMPUTATION  OF  PROPERTIES  FOR 
C     OXYGEN 

C  TO, HO, SO  ARE  FROM  NBS  TN270  (SELECTED  VALUES  OF  CHEM  THERMO  PROPS) 
C  PAGE  11.  COMPUTED  IDEAL  GAS  PROPS  CONFORM  TO  JANAF( NSRDS-NBS37(71) ) 

T0=298.15  $  S0=205.029   $  H0=8680.1 

GI(10)=  HO-HI(TO) 

GI(11)=  SO-SI(TO) 

R=8.20539E-2*. 101325 

DTP=.4081997364372E+02 

DTPV=.3318894767078E-03 

TCC=VP(8) 

PCC=VPN(TCC) 

PTP=VP(9) 

TTP=VP(7) 

TUL=400. 

TLL=TTP 

PUL=120. 

DCC=13.63 

RETURN 

END 
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Appendix  F.  Properties  Subroutines 


SUBROUTINE  PROPS(PP,DD,TT) 
C  THE  32  TERM  EQUATION  OF  STATE,  INPUT  IS  DENSITY(MOLES/L) , 
C  TEMPERATURE(K),  OUTPUT  (PP)  IS  PRESSURE(MPA) ,0R  DP/DD  IN 
C  LITER-MP A/MOLE  OR  DP/DT  MPA/K  OR  S,H,OR  CV  AT  ONE  LIMIT  OF 
C  INTEGRATION 

DIMENSION  X(33) 

DIMENSION  B(33),G(32),VP(9) 

EQUIVALENCE  (B,X) 

C0MM0N/DATA/G,R,GAMMA3VP,DTP,PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

COMMON/1/B 

DATA{ID=1) 

DATA(IZ=1) 

1  CONTINUE 
IF(IZ.LE.O)GO  TO  2 
IZ=0 

2  CONTINUE 
D=DD 
P=PP 
T=TT 

GM=GAMMA 
02=0*0 
D3=D2*D 
D4=D3*D 
D5=D4*D 
D6=D5*D 
D7=D6*D 
D8=D7*D 
D9=D8*D 
D10=D9*D 
D11=D10*D 
D12=D11*D 
D13=D12*D 
TS=SQRT  (T) 
T2=T*T 
T3-T2*T 
T4=T3*T 
T5-T4*T 
F=EXP  (GM*D2) 

GO  TO  (100, 200, 300, 400, 500, 600, 700), K 

ENTRY  PRESS 
C     ENTRY  FOR  PRESSURE,  INPUT  IS  DENSITY 
C     AND  TEMP.  IN  MOL/L  AND  K,  OUTPUT  IS  IN  MPA. 


AND  TEMP. 

K=l 

GO  TO  1 


00  CONTI 

B( 

1) 

B( 

2) 

B( 

3) 

B( 

4) 

B( 

5) 

B( 

6) 

B( 

7) 

B( 

8) 

B( 

9) 

B( 

10) 

B( 

11) 

B( 

12) 

B( 

13) 

B( 

14) 

R( 

15) 

B( 

16) 

B( 

17) 

NUE 

=D2*T 

=D2*TS 

=D2 

=D2/T 

=D2/T2 

=D3*T 

=D3 

=D3/T 

=D3/T2 

=D4*T 

=D4 

=D4/T 

=D5 

=D6/T 

=D6/T2 

=D7/T 

-D8/T 


32 


B(18)=D8/T2 

B(19)=D9/T2 

B(20)=D3*F/T2 

B{21)=D3*F/T3 

B(22)=D5*F/T2 

B(23)-D5*F/T4 

B(24)=D7*F/T2 

B(25)=D7*F/T3 

B(26)=D9*F/T2 

B(27)=D9*F/T4 

B(28)=D11*F/T2 

B(29)=D11*F/T3 

B(30)=D13*F/T2 

B(31)=D13*F/T3 

B(32)=D13*F/T4 

IF(ID.GT.O)GO  TO  102 

B(33)=P-R*D*T 

RETURN 
102  P=0 

M=32 

DO  101  1=1, M 
101  P=P+B(I)*G(I) 

P=P+R*D*T 

PP=P 

RETURN 

ENTRY  DPDD 
C     PARTIAL  OF  PRESSURE  WITH  RESPECT  TO 
C     DENSITY  -  SEE  PRESSURE 
C     ENTRY  FOR  UNITS 

K=2 

GO  TO  1 
200  CONTINUE 

F1=2.00*F*GM*D 

F21=3.000*F*D2 

F22=5.000*F*D4 

F23=7.000*F*D6  +F1*D7 

F24=9.000*F*D8  +F1*D9 


+F1*D3 
+F1*D5 


F25  = 
F26= 
1 
2 
3 
4 
5 
6 
7 
8 
9 

10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 


1.00*F*D10+F1*D11 
3.00*F*D12+F1*D13 
=2.00*D*T 
=2.00*D*TS 


00*D 

00*D/T 

00*D/T2 

00*D2*T 

00*D2 


=3.00*D2/T 


=3.00*D2/T2 

=4.00*D3*T 

=4.Q0*D3 

=4.00*D3/T 

=5,00*n4 

=6.00*D5/T 

=6.00*D5/T2 

=7.00*D6/T 

=8.00*D7/T 

=8.00*D7/T2 

=9.00*D8/T2 

=F21/T2 

=F21/T3 

=F22/T2 

=F22/T4 

=F23/T2 

=F23/T3 

=F24/T2 


33 


B(27)=F24/T4 

B(28)=F25/T2 

P(29)=F25/T3 

B(30)=F26/T2 

B(31)-F26/T3 

B(32)=F26/T4 

M=32 

IF(ID.GT.O)GO  TO  202 

B(33)=P-R*T 

RETURN 

202 

P=0 

DO  201  1=1, M 

201 

P=P+B(I)*G(I) 

P=P+R*T 

PP=P 

RETURN 

ENTRY  DPDT 

C 

PARTIAL  OF  PRESSURE  WITH  RESPECT 

C 

TO  TEMPERATURE  -  SEE  PRESSURE 

C 

ENTRY  FOR  UNITS 

K  =  3 

GO  TO  1 

300 

CONTINUE 

X(  1)=D2 

X(  2)=D2/(2.00*TS) 

X(  3)=0 

X(  4)=-D2/T2 

X(  5)=-2.00*D2/T3 

X(  6)=D3 

X(  7)=0 

X(  8)=-D3/T2 

X(  9)=-2.00*D3/T3 

X(10)=D4 

X(11)=0 

X(12)=-D4/T2 

X(13)=0 

X(14)=-D6/T2 

X(15)=-2.00*D6/T3 

X(16)=-n7/T2 

X(17)=-D8/T2 

X(18)=-2.00*D8/T3 

X(19)=-2.00*D9/T3 

X(20)=-2.00*n3*F/T3 

X(21)=-3.00*D3*F/T4 

X(22)=-2.00*D5*F/T3 

X(23)=-4.00*D5*F/T5 

X(24)=-2.00*D7*F/T3 

X(25)=-3.00*D7*F/T4 

X(26)=-2.00*D9*F/T3 

X(27)=-4.00*D9*F/T5 

X(28)=-2.00*D11*F/T3 

X(29)=-3.00*D11*F/T4 

X(30)=-2.00*D13*F/T3 

X(31)=-3.00*D13*F/T4 

X(32)=-4.00*D13*F/T5 

IF(ID.GT.0)G0  TO  302 

X(33)=PP-R*D 

RETURN 

302 

P=0 

no  301  1=1,32 

301 

P=P+G(I)*X(I) 

PP=P+R*D 

RETURN 

ENTRY  DSDN 

C 

PARTIAL  OF  ENTROPY  WITH 
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+CPOS(T) 


RESPECT  TO  THE  G  COEFFICIENTS 

K=4 

GO  TO  1 

400  CONTINUE 
S=S0-R*L0GF(D*R*T/P0)+(DSDN(D)-DSDN(O))*1000, 
G1=F/(2.00*GM) 
G2=( F*D2-2 .00*G1 ) /( 2 .00*GM) 
G3=(F*D4-4.00*G2)/(2.00*GM) 
G4=(F*D6-6.00*G3)/(2.00*GM) 
G5= ( F*D8-S .00*G4 ) / ( 2 . 0O*GM) 
G6= ( F*D10-10 .00*G5 )/( 2 .00*GM) 

=-D 

=-D/(2.00*TS) 

=0.D0 

=+D/T2 

=2.00*D/T3 

=-D2/2.00 

-O.DO 

=D2/(2.00*T2) 

=D2/T3 
I--D3/3.00 

=0.D0 
l=D3/(3.00*T2) 

=0.D0 
l=D5/(5.00*T2) 

=  2.00*D5/(S.00*T3) 

=D6/(6.00*T2) 

=D7/(7.00*T2) 

=2.00*D7/(7.00*T3) 
i=D8/(4.0n*T3) 

=2.00*G1/T3 
I=3.00*G1/T4 

=2.00*G2/T3 

=4.00*G2/T5 

=2.00*G3/T3 
l=3.00*G3/T4 

=2.nO*G4/T3 

=4.no*G4/T5 

=2.00*G5/T3 

=3.00*G5/T4 

-2.00*G6/T3 

=3.00*G6/T4 

=4.00*G6/T5 
IF(ID.GT.0)G0  TO  402 
RETURN 
402  P=0 

DO  401  1-1,32 

401  P-P+G(I)*X(I) 
PP-P 

RETURN 

ENTRY  DUDN 
C     TERMS  NEEDED  FOR  ENTHALPY  CALCULATION 

K=F 

GO  TO  1 
500  CONTINUE 
C     H=H0+(T*DSDN(D)-DSDN(0))*1000.+(DUDN(D-DUDN(0))*1000.+CP0H(T) 
C     +(P/D-P*T)*1000. 

G1-F/(2.00*GM) 

G2=(F*D2-2.00*G1)/(2.00*GM) 

G3=(F*D4-4.00*G2)/(2.00*GM) 

G4=(F*D6-6.00*G3)/(2.00*GM) 

G5=(F*D8-8.00*G4)/(2.00*GM) 

G6=(F*D10-10.00*G5)/(2.00*GM) 

X(  1)=D*T 

X(  2)=D*TS 


X 

[  1) 

X 

[   2) 

X 

[   3) 

X 

[   4) 

X 

[    5) 

X 

6) 

X 

7) 

X 

8) 

X 

9) 

X 

10) 

X 

11) 

X 

12) 

X 

13) 

X 

14) 

X 

15) 

X( 

16) 

X 

17) 

X 

18) 

X 

19) 

X 

20) 

X 

21) 

X 

22) 

X 

23) 

X 

24) 

X 

25) 

x( 

26) 

X 

27) 

X( 

28) 

X 

29) 

X 

30) 

X 

31) 

X 

32) 
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X(    3 

=D 

Xf    4 

=D/T 

X(    5 

-D/T2 

X{    6 

)=D2*T/2.00 

X(    7 

)=D2/2.00 

X(    8 

)=D2/(2.00*T) 

X(    9 

)=D2/(2.00*T2) 

X(10 

)=D3*T/3.00 

X(ll 

=D3/3.00 

X(12 

=D3/(3.00*T) 

X(13 

)=D4/4.00 

X(14 

=D5/(5.00*T) 

X(15 

=D5/(5.00*T2) 

X(16 

)=D6/(6.00*T) 

X(17 

=D7/(7.00*T) 

X(18 

)=D7/(7.00*T2) 

X(19 

)=D8/(8.00*T2) 

X(20 

=G1/T2 

X(21 

=G1/T3 

X(22 

=G2/T2 

X(23 

=G2/T4 

X(24 

=G3/T2 

X(25 

-G3/T3 

X(26 

=G4/T2 

X(27 

=G4/T4 

X(28 

=G5/T2 

X(29 

)=G5/T3 

X(30 

=G6/T2 

X(31 

=G6/T3 

X(32 
IF(I[ 

=G6/T4 
).GT.0)G0  TO  502 

RETU? 

m 

502    P=0 

DO   5( 

n  1=1,32 

501    P=P+C 

un*x(i) 

PP  =  P 

RETUf 

ni 

EMTR^ 

(   TDSDT 

C            TEMP 

TIMES   THE    PARTIAL 

OF 

C             EMTR( 

)PY  WITH   RESPECT 

TO 

TEMP. 

K=6 

GO   T( 

)    1 

600  COMT 

[HUE 

c        cv=c\ 

/0+(TDSDN(/)-TDSDM(D))*1000 

G1=F/ 

a2.00*GM) 

G2=(F 

■*D2-2.00*Gl)/(2 

00- 

^GM) 

G3  =  (F 

-*D4-4.00*G2)/(2 

00^ 

^GM) 

G4=(F 

-*D6-6.00*G3)/(2 

00^ 

^GM) 

G5-(F 

■*D8-8.00*G4)/(2 

00*GM) 

G6=(F 

-*Din-10.00*G5)/(2.( 

30*GM) 

X(l)  = 

=0 

X(    2 

=-D/(4.00*TS) 

X(3)^ 

=0 

X(    4 

=2.00*D/T2 

X(    5 

=6.00*D/T3 

X(6)  = 

=0 

X(7)=0 

X(    8)=D2/T2 

X(    9)=3.00*D2/T3 

X(10)=0 

X(11)=0 

X(12)=(2.00*D3)/(3.00*T2) 

X(13)=0 

X(14)=(2.0n*D5)/(5.00*T2) 

X(15)=(6.00*D5)/(5.00*T3) 
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X(16)=D6/(3.00*T2) 

X(17)-(2.00*D7)/(7.00*T2) 

X(18)=(6.00*D7)/(7.00*T3) 

X(19)-(3. 00*08)7(4. 00*T3) 

X(20)=6.000*G1/T3 

X(21)=12.00*G1/T4 

X(22)=6.000*G2/T3 

X(23)=20.00*G2/T5 

X(24)=6.000*G3/T3 

X(25)=12.00*G3/T4 

X(26)-6.00n*G4/T3 

X(27)=20.00*G4/T5 

X(28)-6.000*G5/T3 

X(29)=12.00*G5/T4 

X(30)=6.00n*G6/T3 

X(31)=12.00*G6/T4 

X(32)=2n.00*G6/T5 

IF(ID.GT.0)G0  TO  602 

RETURN 
602    P=0 

DO  601    1=1,32 
601    P=P+G(I)*X(I) 

PP=P 

RETURN 

ENTRY   DP2D2 

SECOND    PARTIAL    OF    PRESSURE    WITH 

RESPECT  TO   DENSITY    SQUARED 

K=7 

GO   TO    1 
700  CONTINUE 

F1=2.*F*GM*D 

F12=2.*F1*GM*D+2.*F*GM 

F212=3.*F1*D2+3.*2.*D*F+F12*D3+F1*3.* 

F222=5 .*F1*D4  +5 .*4 .*D3*F+5 .*D4*F1+F1 

F232=7  .*Fl*D6+7 .*6 .*D5*F+7 .*D6*F1+F12 

F242=9 .*F1*D8+P .*8 .*D7*F+9 .*D8*F1+F12 

F252=11.*F1*D10+10.*11.*D9*F+11.*D10* 

F262=13.*F1*D12+13.*12.*D11*F+13.*D12 

B(1)=2.*T  $B(2)=2.*TS  $  B(3)=2. 

B(4)=2./T  $  B(5)=2./T2  $  B(6)=6.* 

B(7)=6.*D  $  B(8)=6.*D/T  S  B(9)=6 

B(10)=12.*D2*T  $  B(11)=12.*D2  $  B 
B(14)=30.*D4/T  S  P 
S  B(17)=56.*D6/T  $ 
%     B(20)=F212/T2  S 


702 


701 


B(13)=20.*D3  $ 

B(16)=42.*D5/T 

B(19)=72.*D7/T2 

B(22)=F222/T2 

B(23)=F222/T4  $ 

R(26)=F242/T2  $ 

B(29)=F252/T3  $ 

B(32)=F262/T4 

M=32 

IF(ID.GT.O)GO  TO  702 

B(33)=PP 

RETURN 

P=0 

DO  701  1=1, M 

P=P+B(I)*G(I) 

PP=P 

RETURN 

END 


B(24)=F232/T2 
B(27)=F242/T4 
B(30)=F262/T2 


D2 

2*05 

*07 

*D9 

F1+F12*D11 

*F1+F12*D13 

0*T 

*D/T2 

(12)=12.*02/T 

(15)=30.*D4/T2 
B(18)=56.*D6/T2 
B(21)=F212/T3 

$  B(25)=F232/T3 
$  B(28)=F252/T2 
$  B(31)=F262/T3 


FUNCTION  VPN(TT) 
CALCULATES  VAPOR  PRESSURE(MPA] 


INPUT  IS  TEMP(K) 
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DIMENSION  G(32),VP(9) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 
T=TT 

X=(1.-VP(7)/T)/(1.-VP(7)/VP(8)) 

VPN=VP(9)*EXP  (VP(1)*X+VP{2)*X*X+VP(3)*X**3+VP(4)*X**4+VP(5)*X- 
1(1.-X)**VP(6)) 
RETURN 
END 


FUNCTION  FINDTV(POBS) 
C   ITERATES  VAPOR  PRESS  EON  TO  FIND  TEMP(K),  FOR  INPUT  OF  PRESS(MPA). 
C     GIVEN  AN  INPUT  PRESSURE( MPA) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

DIMENSION  G(32),VP(9) 

T=VP{8) 

DO  7  1=1,10 

P=VPN(T) 

IF(ABS  (P-P0BS)-.000001*P0BS)8,8,6 

6  CONTINUE 
CORR=(PnBS-P)/DPDTVP(T) 

7  T=T+CORR 

8  CONTINUE 
FINDTV=T 
RETURN 
END 


FUNCTION  CV(D,T) 
CALCULATES  CV( J/( MOL*K) ) .   INPUT  DENS{MOL/L)  AND  TEMP(K). 
DATA(R=8. 31434) 
DD=D 
TT=T 

CALL.  TDSDT(CD,DD,TT) 
DD=0 

CALL  TDSDT(C0,DD,TT) 
CV=CPI(TT)+(C0-CD)*1000. 
CV=CV-R 
RETURN 
END 


FUNCTION  FIND  D(P,T) 
C  ITERATES  EOUATION  OF  STATE  FOR  DENSITY( MOL/L) ,  FOR  GIVEN  PRESSURE(MPA)  AND 
C  TEMPERATllRE(K) .   IF  ITERATION  FAILS  TRY  USING  FUNCTION  FIND  m, 

DIMENSION  G(32),VP{9) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

TT=T 

IF(TT.GT.VP(8)*.99999)G0  TO  100 

IF(  P.GT.VPN(TT))G0  TO  101 

DD=SATV(TT) 


GO   TO    102 

100 

PC=PCC 

X=(1.1/(9.*PC)) 

DD=P/(R*T*X) 

*P+ 

.7/9. 

IF(P/PC.GT.20.. 

AND 

.T/VP(8).LT 

.2 

5)DD= 

=DTP 

GO   TO   102 

101 

DD=SATL(TT) 

102 

CONTINUE 

DO   10    1=1,50 

IF(DD.LE.0.0.0R 

.DD 

.GT.50.)G0 

TO 

11 

CALL  PRESS(PP,DD,TT) 
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IF(PP.LE.O.O)GO   TO   11 
P2=PP 

IF(ABS   (P-P2)-1.E-7*P)20,20,1 
1    CALL   DPDD(PP,DD,TT) 
DP-PP 

C0RR=(P2-P)/DP 
IF(ABS   (C0RR)-1.E-7*DD)20,20,10 

10  DD=DD-CORR 

11  CALL    REGULA(P,DD,T) 
20   FIND    D=DD 

RETURN 
END 


SUBROUTINE    REGULA( PI ,DD ,TT) 
ITERATES   EON   OF   STATE   FOR   DENSITY   WHEN   SUBPROG   FINDD   FAILS. 
DIMENSION  G(32)  ,VP(9) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 
T=TT 
P=PI 
D2-0 

IF(T.LT.TCC)GO   TO   10 
DO=DCC*TCC/T 
GO  TO   20 

10  PP=VPN(T) 
IF(P.GT.PP)GO   TO   15 
DO-SATV(T) 

DO   11    1-1,150 

CALL    PRESS(P0,D0,T) 

IF(P0.GE.P)G0   TO    12 

11  D0-D0+.0001*DO 
GO   TO  42 

12  D1=D0 

13  CALL   PRESS(P1,D1,T) 
IF(P1.LT.P)G0   TO   14 
IF(D1.LE..1*PTP)G0   TO   42 
D0=D1 

Z=(P1-P)/P 
IF(Z.LT..l)Z=.l 
IF(Z.GT..9)Z=.9 
D1=D1-Z*D1 
GO  TO   13 

14  CALL    PRESS(PO,DO,T) 
DO   140    1=1,50 

D=D1 

P3=P1 

IF(ARS(P-P1) .LT..00001*P)G0   TO   40 

P2-P-P1 

D1=D1+(D1-D0)*P2/(P1-P0) 

IF(ABS(D-Dl).LE..nn001*D)G0   TO  40 

IF(ABS(P-P1).LT..005*P)D2=FIND   M(P,T,D1) 

IF(D2.GT.O.O.AND.02.LT.50.)D1=D2 

D2=0 

CALL   PRESS(P1,D1,T) 

IF{P0.GT.P.AND.P1.GT.P)G0   TO    120 

IF(P0.LT.P.AND.P1.LT.P)G0   TO    120 

GO   TO    140 
120   P0=P3 

DO=D 
140  CONTINUE 

GO  TO  41 

15  DO=SATL(T) 
DO   16    1=1,10 

CALL    PRESS(PO,DO,T) 
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IF(PO.LE.P)GO  TO   17 

16  D0=D0-.0001*D0 
GO  TO  42 

17  Dl-DO 

18  CALL    PRESS(P1,D1,T) 
IF(D1.GE.50.)G0  TO  42 
IF(P1.GT.P)G0  TO   14 
D0=D1 

Z=(P-P1)/P 

Z=Z*10 

IF(T/TCC.LT. .6)7=1. 

IF(Z.LT.1.)Z=1. 

IF(Z.GT.9.)Z=9. 

D1=D1+.01*D1*Z 

GO  TO   18 

20  CALL  PRESS(Pn,DO,T) 
IF(P.LE.PO)GO  TO  30 
D1=D0 

21  CALL   PRESS(P1,D1,T) 
IF(P1.GE.P)G0   TO   14 
IF(D1.GE.50.)G0   TO   42 
D0=D1 

Z=(P-P1)/P 

Z=Z*10 

IF(Z.LT.1)Z=1 

IF(Z.GT.9)Z=9 

Dl=ni+.1*D1*Z 

GO  TO   21 

30  D1=D0 

31  CALL    PRESS(P1,D1,T) 
IF(P1.LE.P)G0   TO    14 
IF(D1.LE..1*PTP)G0   TO   42 
D0=D1 

Z=(P1-P)/P 
Z=Z*10 

IF{Z.LT.1)Z=1 
IF(Z.GT.9)Z=9 
D1=D1-.1*D1*Z 
GO  TO   31 
40   DD=D1 


RETURN 
41  PRINT  101,P,T,D 
102  FORMAT{*  REGULA  FAILED  AT  P=*,F7.2,*  AND  T=*,F7.2) 
101  F0RMAT(*  DENSITY  ITTERATION  FAILED  AT  P=*,F7.2,*  AND  T=*,F7.2, 
1/*  DENSITY  RETURNED  IS*,E17.8) 
RETURN 


42  PRINT  102, P,T 
RETURN 
END 


FUNCTION  CP(D,T) 
CALCULATES  CP( J/( M0L*K ) ) .   INPUT  DENS(M0L/L),  TEMP(K) 
CVEE=CV(D,T) 
CALL  DPDT(DPT,D,T) 
CALL  DPDD(DPD,D,T) 

CP=CVEE+(T/(D**2)*(DPT**2)/DPD)*1000. 
RETURN 
END 
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FUNCTION  DPDTVP(TT) 
C  CALCULATES  THE  DERIVATIVE  OF  PRESSURE  WITH  RESPECT  TO  TEMPERATURE 
C  AT  SATURATION.   INPUT  IS  TEMP(K),  OUTPUT  IS  DPDT(MPA/K). 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

DIMENSION  G(32),VP(9) 

T=TT 

IF(TT.GT.VP(8))G0  TO  1 

X=(1.-VP(7)/T)/(1.-VP(7)/VP(8)) 

DXDT=(VP(7)/T**2)/(1.-VP(7)/VP(8)) 

DPDT=VP(1)*DXDT+2.*VP(2)*X*DXDT+VP(3)*3.*X**2*DXDT+VP(5)* 
1((1.-X)**VP(6))*DXDT+VP(5)*X*((1.-X)**(VP(6)-1.))*VP(6)*(-DXDT) 

DPDT=DPDT*VPN(T) 

DPDTVP=DPDT 

RETURN 
1  DPDTVP=0 

RETURN 

END 


FUNCTION  FIND  M{P,T,DD) 
C  SOLVES  FOR  DENSITY(MOL/L)  BY  ITERATION.   INPUT  IS  PRESSURE(MPA) , 
C  TEMPERATURE(K),  AND  A  STARTING  VALUE  OF  DENSITY.   THIS  FCN  IS  AN 
C  ALTERNATIVE  FOR  FUNCTION  FIND  D. 

TT=T 

DO  10  1=1,50 

CALL  PRESS(PP,DD,TT) 

P2=PP 

IF(ABS  (P-P2)-1.E-7*P)20,20,1 
1  CALL  DPDD(PP,DD,TT) 

DP=PP 

C0RR=(P2-P)/DP 

D=DD 

IF(ABS  (C0RR)-1.E-7*D)20,20,10 
10  DD=DO-CORR 

FIND  M=0 

RETURN 
20  FIND  M=DD 

RETURN 

END 


FUNCTION  ENTHAL(P,D,T) 
C  RETURNS  ENTHALPY(J/MOL)   FOR  AN  INPUT  OF  PRESSURE(MPA)  AND  DENSITY(MOL/L] 
C  AND  TEMPERATURE(K). 

R=  .00831434 

DD=D 

TT=T 

CALL  DSDN{SD,DD,TT) 

CALL  DUDN(UD,DD,TT) 

DD=0 

CALL  DSDN(SO,DD,TT) 

CALL  DUDN(UO,DD,TT) 

ENTHAL=T*(SD-S0)*1000.+(UD-U0)*1000.+HI{T)+(P/D-R*T)*1000. 

RETURN 

END 


FUNCTION  ENTROP(D,T) 
CALCULATES  ENTROPY(J/(MOL-K) ,  FROM  INPUT  OF  DENSITY(MOL/L)  AND  TEMP(K) 
R=  .00831434 
P0=  .101325 
DD=D 
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TT=T 

CALL  DSDN(SD,DD,TT) 

DD=0 

CALL  DSDN(SO,DD,TT) 

ENTROP=(SD-S0)*1000.-P,*ALOG(D*R*T/P0)*1000.+SI(T) 

RETURN 

END 


FUNCTION  SATL(TT) 
CALCULATES  DENSITY(MOL/L)  OF  SATURATED  LIQUID.   INPUT  IS  TEMP(K) 
DIMENSION  A(20) 
DIMENSION  G(32),VP(9) 

C0MM0N/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 
COMMON/SATC/A,DTPV 
T=TT 
K  =  14 
KK=7 

GO   TO   10 
ENTRY   SATV 
K  =  l 
KK=13 
T=TT 

10  IF(T.GE.TCC*.995^5?)G0   TO   20 
ITT=TCC 

IF(ITT+1-T.LT.1.)T=ITT 
X-(T-TCC)/(TTP-TCC) 
D=A(K)*ALOG(X) 

DO   11    1=2, KK 
K=K  +  1 
MM=I 
IF(MM.GE.5)MM=MM+1 

11  D=D+A(K)*(l.-X**((MH-5)/3.)) 
IF(K.LT.14)G0  TO   12 
D=DCC+EXP(D)*(DTP-DCC) 

GO  TO   13 

12  D=DCC+EXP(D)*(DTPV-DCC) 

13  SATL=D 
IF(ITT+1-TT.LT.1.)SATL=D-(D-DCC)*(TT-T) 
RETURN 

20  DSATL=DCC 
RETURN 
END 


FUNCTION  SOUND(D,T) 
CALCULATES  SPEED  OF  SOUND(M/S).   INPUT  IS  DENSITY(MOL/L)  AND  TEMP(K). 
COMMON/CRIT/  W,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 
CALL  DPDD(DP,D,T) 

S0UND=((CP(D,T)/CV(D,T))*DP*1000000./W)**.5 
RETURN 
END 


FUNCTION  VISC(DD,T) 
CALCULATES  VISCOSITYf  MICRO  PA*S)  .   INPUT  IS  DENSITY(MOL/L)  AND  TEMP(I<) 
COMMON/CRIT/  EM,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 
D=DD*EM/1000. 

VISC=DILV(T)+FDCV(D,T)+EXCESV(D,T) 
RETURN 
END 
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FUNCTION  THERM(DD,T) 
RETURNS  THERMAL  CONDUCTIVITY( W/( M*K) ) .   INPUT  IS  DENSITY(MOL/L) ,  TEMP(K) 
COMMON/HAN/CR,TCI 
COMMON/ I SP/N,NW,NWW 

COMMON/CRIT/  EM,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 
D=DD*EM/1000. 
CR=CRITC(n  T) 

THER=DILT(f)+FDCT(D,T)+EXCEST(D,T)+CR 
TCI=THER-CR 
THERM=THER 
RETURN 
END 


C 

FUNCTION  EXCESV(D,T) 
C     CALCULATES  EXCESS  VISCOSITY 

COMMON/DATAl/GV ,GT ,FV , FT , EV , ET 

DIMENSION  GV(9),GT(9),FV(4),FT(4),EV(8),ET(8) 

R2=D**(.5)*((D-EV(8))/EV(8)) 

R=D**(.l) 

X=EV(1)+EV(2)*R2+EV(3)*P+EV(4)*R2/(T*T)+EV(5)*R/T**(1.5)+EV(6)/T 
1+EV(7)*P,2/T 

X1=EV(1)+EV(6)/T 

EXCESV=  EXP(X)-  EXP(Xl) 

RETURN 

ENTRY  EXCEST 
C     CALCULATES  EXCESS  THERMAL  CONDUCTIVITY 

R2=D**(.5)*((D-ET(8))/ET(8)) 

R=D**(  .1) 

X-ET(1)+ET(2)*R2+ET(3)*R+ET(4)*R2/(T*T)+ET(5)*R/T**(1.5)+ET(6)/T 
1+ET(7)*R2/T 

X1=ET(1)+ET(6)/T 

EXCESV=  EXP(X)-  EXP(Xl) 

RETURN 

END 


FUNCTION  FDCV(D,T) 
FIRST  DENSITY  CORRECTION  FOR  VISC  AND  THERMAL  COND. 
COMMON/DATAl/GV, GT,FV, FT, EV,ET 
DIMENSION  GV(9) ,GT(9) ,FV(4) ,FT(4)  ,EV(8)  ,ET(8) 
FDCV=(FV(1)+FV(2)*(FV(3)-AL0G(T/FV(4)))**2)*D 
RETURN 
ENTRY  FDCT 

FDCV=(FT(1)+FT(2)*(FT(3)-AL0G(T/FT(4)))**2)*D 
RETURN 
END 


FUNCTION  CRITC(D,T) 
C  CALCULATES  CRITICAL  ENHANCEMENT  FOR  THERM.  COND. 
C  INPUT  UNITS  ARE  G/CC,  K,  OUTPUT  IS  W/(M*K). 

COMMON/CRIT/  EM,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 

COMMON/CHECK/DELD,DELT,DSTAR,TSTAR 

COMMON/HJM/EPSI,CPCV,RRR,AKT 

COMMON/ I SP/N,NW,NWW 

AV=6.0225E+23  $  BK=1 .38054E-16 

DELD=ABS  (D-DC)/DC  $  DELT=ABS  (T-TC)/TC 
C  CALCULATE  DISTANCE  PARAMETER 

R=(RM**2.5)*(D**O.5)*(AV/EM)**0.5 

R=R*(E0K**0.'^)*X/(T**0.5) 
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RRR=R 
C   GENERAL  EQUATION 

DX=D*1000.0/EM 
C   DX  IN  MOL/L,   D  IN  G/CM3. 

CALL  DPDT(DPT,DX,T) 
C  DPDT  IN  MPA/K. 

DPT=DPT*1.0E+7 
C   DPDT  NOW  IN  DYNES/{CM2*K) 

CALL  DPDD(DPD,DX,T) 
C   DPDD  IN  L*MPA/MOL. 

DPD=DPD*1.0E+7*1000./EM 
C  DPDD  NOW  IN  DYNE*CM/G. 

IF(  DPD.LT.0.0)  DPD^l.O 
94  VIS=VISC(DX,T)*(1.0E-05) 
C  VISCOSITY  NOW  IS  G/{CM*S). 

IF(DELD.E0.0.25.  OR.  DELD.LT.  0.25)  8,10 

8  IFCDELT.EQ. 0.025.    OR.    DELT.LT. 0.025    )    9,    10 

9  COMPRES=SEMG{D,T) 
GO   TO    12 

10  C0MPRES=1.0/(D*DPD)**0.5 
12   EX=BK*T**2*(DPT**2)*C0MPRES 

EXB=R*((BK*T)**0.5)*(D**0.5)*((AV/EM)**0.5) 

CRIT=EX/(EXB*6.0*3.14159*VIS) 
C   THERMAL  COND,  CRIT,  IS  IN  ERG/( CM*SEC*K) 
C   PUT  IN  DAMPING  FACTOR 

BDD=((D-DC)/DC)**4 

BTT-((T-TC)/TC)**2 

BXX-  -18.66*PTT  -  4.25*PDD 

IF(BXX.LT.-l.E+2)  BXX=  -l.E+2 

FACT=  EXP(  BXX  ) 
C     FACT=EXP  (-18.66*BTT  -  4.25*BDD) 

DELC=CPIT*FACT 

CRITC=DELC/100000. 
C   THERMAL  COMD,  CRITC,  IS  NOW  IN  W/(M*K) 

AKT=COMPRES*COMPRES 

EPSI  =  R*T^*BK*T*(AV*D/EM)*AKT 

EPSI=EPSI**0.5 
C   CALC  CP-CV 

CPCV=T*(DPT**2)*AKT/D 

RETURN 

END 


FUNCTION  SENG(D,T) 
SCALED  EOUATION  OF  STATE  FOR  CRITICAL  REGION 
COMMON/CRIT/  EM,  EOK,  RM,  TC,  DC,  X  ,  PC,  SIG 
COMMON/SEN/BETA, XO, DELTA, El,  E2,  AGAM 
COMMON/CHECK/DELD,DELT,DSTAR,TSTAR 
DSTAR=  D/DC  $  TSTAR-T/TC 
BET0=1./BETA 
XX=DELT/DELD**RETO 
AG=AGAM-1.0 

BET2=  2.0*BETA 
AGP=AG/BET2 
DEL1=DELTA-1.0 

AGBB=(AG-BET2)/BET2 
XXO=(XX+  XO)/XO 
XXB=XX0**BET2 
BRAK=1.0  +  E2*XXB 

BRAK1=BRAK**AGB 
H=E1*XX0*BRAK1 

HPRIM=(E1/X0)*BRAK1  +  ( AG/X0)*E1*E2*(XXB)*( BRAK**AGBB) 
RC0M=(DELD**DEL1)*(DELTA*H  -  (XX/BETA)*HPRIM  ) 
RCOMP= 1 . 0/ ( RC0M*DSTAR**2 ) 
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c 


RCM=RC0MP/(PC*1.0E+7) 
RCM    IN   CM2/DYNIE,    PC    IN   MPA 
RCM=RCM**0.5 

.SENG=RCM 
RETURN 
END 


FUNCTION  DILV(T) 
DILUTE  GAS  VISCOSITY  AND  THERMAL  CONDUCTIVITY. 

C0MM0N/DATA1/GV,GT,FV,FT,EV,ET 

DIMENSION  GV(9),GT(9),FV(4),FT(4),EV(8),ET(8) 

SUM=0 

TF=T**(l./3.) 

TFF-T**(-4./3.) 

DO  10  1=1,9 

TFF=TFF*TF 
10  SUM=SUM+GV(I)*TFF 

DILV=SUM 

RETURN 

ENTRY  DILT 

TF=T**(l./3.) 

TFF=T**(-4./3.) 

SUM=0 

DO  20  1=1,9 

TFF=TFF*TF 
20  SUM=SUM+CT(I)*TFF 

DILV=SUM 

RETURN 

END 


FUNCTION  FIND  P(D,T) 

DIMENSION  G(32)  ,VP(9) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

DD=D 

TT=T 

IF(TT.LT.TCC)GO  TO  10 
1  CALL  PRESS(PP,DD,TT) 

FIND  P=PP 

RETURN 
10  P=VPN(TT) 

DV=FIND  D(P-.n001,TT) 

DL=FIND  D(P+.0001,TT) 

IF(DD.LE.DV.OR.DD.GE.DL)G0  TO  1 

PRINT  100,DV,DL,DD 

CALL  PRESS(PP,DV,TT) 

FIND  P=PP 

D=DV 

RETURN 
100  FORMAT(*  THE  STATE  POINT  YOU  HAVE  SPECIFIED  CORRESPONDS  TO  A  * 
1/*  DENSITY  IN  THE  LIQUID  VAPOR  COEXISTENCE  REGION* 
2/*  THE  DENSITY  OF  THE  SATURATED  VAPOR  IS  *,F6.4,*  MOLES/LITER* 
3/*  THE  DENSITY  OF  THE  SATURATED  LIOUID  IS  *,F8.4,*  MOLES/LITER* 
4/*  AND  THE  INPUT  DENSITY  IS  *,F8.4,*  MOLES/LITER* 
5/*  SATURATED  VAPOR  IS  ASSUMED*) 

END 


FUNCTION  FINDT(P,D) 
RETURNS  TEMPERATURE(K)  ,  FROM  THE  32-TERM  MBWR  EON  OF  STATE, 


45 


INPUT  IS  PRESSURE(MPA)  AND  DENSITY( MOL/L) . 
DIMEMSIOM  G(32),VP(9) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 
PP=P 
DD=D 

IF(P.GE.PCC)GO  TO  1 
TSAT=FINDTV(PP) 
DV=FIMD  D(PP-. 00001, TSAT) 
DL=FIMD  D(PP+. 0001, TSAT) 
IF{nD.GT.DV.AND.DD.LT.DL)GO  TO  30 
TT=TSAT 
GO  TO  2 

1  TT=TCC 

2  DO  10  1=1,10 

CALL  PRESS(P2,DD,TT) 
IF(ABS(PP-P2)-1.E-7*PP)20,20,11 
11  CALL  nPDT(DP,DD,TT) 
C0RR=(P2-PP)/DP 
IF(ARS(C0RR)-1.E-'^)20,20,10 
10  TT=TT-CORR 
20  FINDT=TT 

RETURN 
30  FINDT=TSAT 
D=DV 

PRINT  100,DV,DL,DD 
100  F0RMAT(*  THE  STATE  POINT  YOU  HAVE  SPECIFIED  CORRESPONDS  TO* 
1/*  A  DENSITY  IN  THE  LIQUID  VAPOR  COEXISTENCE  REGION* 
2/*  DENSITY  OF  THE  SATURATED  VAPOR  IS*,F8.4,*  MOLES/LITER* 
3/*  DENSITY  OF  THE  SATURATED  LIOUID  IS*,F8.4,*  MOLES/LITER* 
4/*  INPUT  DENSITY  IS*,F8.4,*  MOLES/LITER* 
5/*  SATURATED  VAPOR  CONDITIONS  ARE  ASSUMED*) 
RETURN 
END 


SUBROUTINE  INFO 

PRINT  100 
100  F0RMAT(/*  WHEN  THE  PROGRAM  ASKS  FOR  A  FLUID  SELECTION,  ENTER  THE* 
1*APPR0PRIATE  NUMBER.*/*  AN  INAPPROPRIATE  NUMBER  WILL  TERMINATE* 
2*  THE  PROGRAM.  WHEN  THE  PROGRAM  ASKS*/*  FOR  A  PRESSURE,  DENSITY* 
3*,  AND  TEMPERATURE,  ENTER  ANY  TWO  OF  THE  THREE*/*  AND  A  ZERO  FOR* 
4*  THE  THIRD.*/*  THE  ORDER  MUST  BE  P,D,T,  AND  ONE  OF  THE  THREE* 
5*  MUST  BE  ZERO.*/*   IF  ALL  THREE  ARE  ZERO  THE  PROGRAM  ASKS  FOR  A* 
6*  NEW  FLUID.*) 

RETURN 

END 


FUNCTION  DIEL(P,D,T) 
DIELECTRIC  CONSTANT.   INPUT  P(MPA),  D(MOL/L)  AND  T(K). 
C0MM0N/DIEL/BX(6)  ,PX(6) 

CM=  BX(1)+  BX(2)*D+  BX(3)*D**2+  BX(4)*D**3+  BX(5)*P+  BX(6)*T 
DIEL=(1.+2.*D*CM)/(1.-D*CM) 
RETURN 
END 


FUNCTION  CPI(T) 
C   CALCULATES  SPECIFIC  HEAT,  ENTROPY,  AND  ENTHALPY  FOR  THE  IDEAL  GAS. 
C  OUTPUT  IS  IN  J/(MOL*K),  FOR  CP  AND  S,   AND  J/MOL  FOR  H. 
C  HYDROGEN,  (N=l),  IS  TREATED  AS  A  SPECIAL  CASE  AS  THE  COEFF.  FOR 
C  CP  ARE  IN  THREE  TEMPERATURE  RANGES.   T  <  40  K ,   40  <  T  <  140  K, 
C  AND  T  >  140  K. 
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COmON/CPID/G(ll)  ,GH(11)  ,GL(11) 

COMMON/ISP/N,MW,NWW 

K  =  l 

IF(H.EQ.l)  10,15 
10  C=CPO(T) 
15  U=G(9)/T 

EU=EXP  (U) 

TS=1./T**4 

GO  TO  (20,40,55),K 
20  CPI-G(8)*U*U*EU/(EU-1.)**2 

DO  25  1=1,7 

TS=TS*T 
25  CPI=CPI+G(I)*TS 

CPI  =  CPI*8. 31434 

RETURN 

ENTRY  SI 

K=2 

IF(N.EO.l)  30,35 
30  C=CPO(T) 
35  GO  TO  15 
40  CPI=G(8)*(U/(EIJ-1.)-AL0G(1.-1./EU)) 

l-G(l)*TSn/3.-G(2)*TSn*T/2.-G(3)/T+G(4)*AL0G(T)+G(5)*T+G(6)*Tn/2 
2.+G(7)*T**3/3. 

CPI=CPI*8.31434+G(11) 

RETURN 

ENTRY  HI 
K=3 

IF(N.EQ.l)  45,50 
45  C=CPO(T) 
50  GO  TO  15 

55  CPI=G(8)*Un/(EU-l.)-G(l)/(2.*T*T)-G(2)/T+G(3)*AL0G(T)+G(4)*T 
l+G(5)*T*T/2.+G(6)*T**3/3.+G(7)*T**4/4. 

CPI  =  CPI*8. 31434+0(10) 

RETURN 

END 


FUNCTION  CPO(T) 
C  SELECTS  PROPER  COEFF  FOR  G(I),  AND  COMPUTES  G(10)  AND  G(ll) 
C  FOR  HYDROGEN  ONLY. 

C0MM0N/CPID/G(11) ,GH(11) ,GL(11) 

TX1=140. 

TX2=40. 

DO  10  J=l,ll 
10  G(J)=GH(J) 

IF(T.LT.140.)  30,20 
20  CP0=1.0 

RETURN 
30  G(10)=G(10)+GHI(TX1) 

G(11)=G(11)+GSI(TX1) 

DO  40  J=l,8 
40  G{J)=GL(J) 

G(10)=G(10)-GHI(TX1) 

G(11)=G(11)-GSI(TX1) 

IF(T.LT.40.)  60,50 
50  CP0=1.0 

RETURN 
60  G(10)=G(10)+GHI(TX2) 

G(11)=G(11)+GSI(TX2) 

DO  70  J=l,8 
70  G(J)=0.0 

G(4)=2. 5000315 

G(10)=G(10)-GHI(TX2) 


47 


G(11)-G(11)-GSI(TX?) 

CP0=1.0 

RETURN 

END 


FUNCTION  GHI(T) 

COMMON/CPID/Gdl)  ,GH(11)  ,GL(11) 
1  U=G(9)/T 
EU=EXP(U) 

GHI=G(R)*U*T/(EU-1.)-G(1)/(2.*T*T)-G(2)/T+G(3)*AL0G(T)+G(4)*T 
A  +G(5)*Tn/2.+G(6)*T**3/3.+G(7)n**4/4. 
GHI=GHI*8. 31434 
RETURN 

ENTRY  GSI 

U=G(9)/T 

EU=EXP(U) 

TS=1./T**4 

GHI=  G(8)*(U/(EU-1.)-AL0G(1.-1./EU))- 
A  G(l)*TS*T/3.-G(2)*TS*T*T/2.-G(3)/T+G(4)*AL0G(T)+G(5)*T+ 
B  G(6)n*T/2.+G(7)*T**3/3. 

GHI=GHI*8. 31434 

RETURN 

END 


SUBROUTINE  LIMITS(P,T,IL) 

DIMENSION  G(32),VP{9) 

COMMON/DATA/G,R, GAMMA, VP, DTP, PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC 

C0MM0N/DIEL/BX(6),PX(6) 

IF(P.GT.PUL)GO  TO  10 

IF(T.GT.TUL.OR.T.LT.TLL)GO  TO  12 

PM=PMELT(T)  ' 

IF(P.GT.PM)  GO  TO  20 

IL=1  S  RETURN 

10  PRINT  11,PUL 

11  FORMAT! *  THE  INPUT  PRESSURE  IS  OUT  OF  THE  RANGE  OF  THIS  EQUATION  * 
1/*  THE  RANGE  FOR  THIS  EQUATION  IS  FROM  0  TO  *,F6.0,*  BAR*) 

IL=0  S  RETURN 

12  TLLF=  (TLL-273.15)*1.8+32. 
TULF=  (TnL-273.15)*1.8+32. 
PRINT  13,  TLL,TUL,TLLF,TULF 

13  F0RMAT(*  THE  INPUT  TEMPERATURE  IS  OUT  OF  RANGE* 

A  /*  THE  RANGE  FOR  THIS  EQUATION  IS  *,F6.2,*  K   TO  *,F6.0,*  K*,/, 
B  27X,*  OR  *,F8.2,*  F   TO  *,F6.0,*  F*) 

IL=C  $  RETURN 
20  TM=TMELT(P) 

TF=(TM-273.15)*1.8+32. 

PRINT  22,  TM,TF 
22  FORMAT(*  SOLID  PHASE  DETECTED.*,/,*  FOR  THIS  PRESSURE,  TEMP* 
A  *  SHOULD  EXCEED  *,F8.3,*  K,  0R*,F9.3,*  F*) 

IL=0  $  RETURN 

END 


FUNCTION  PMELT(T) 
COMPUTES  MELTING  PRESSURE(MPA)  FOR  INPUT  TEMPERATURE(K) 

C0MM0N/DIEL/BX(6),PX(6) 

COMMON/ I SP/N,NW,NWW 

IF(N.EO.l)  20,10 
10  PMELT=  PX(1)+  PX(2)*T**PX(3) 

RETURN 
20  IF(T.GT.22.)  30,10 
30  PMELT=  PX(4)+  PX( 5)*T**PX(6) 
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RETURN 
END 


FUNCTION  TMELT(P) 
COMPUTES  MELTING  TEMPERATURE(K)  FOR  INPUT  PRESSURE(MPA) 

COMM0N/DIEL/BX(6),PX(6) 

COMMON/ I SP/N,NW,NWW 

IF(N.EO.l)  20,10 
10  TMELT=((P-PX(1))/PX(2))**(1./PX(3)) 

RETURN 
20  IF(P.GT.31.64)  30,10 
30  TMELT=((P-PX(4))/PX(5))**(1./PX(6)) 

RETURN 

END 
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